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SUMMARY: The responses of systems, whose equations of motion are linear, to 
arbitrary applied forces are conveniently calculable by the use of impulsive 
or indicial admittances. Expressions for these admittances are here given 
explicity for some simple dynamical chains and the procedure for their 
calculation in chains with more degrees of freedom is outlined. Some general 
characteristics of the dynamical chains are also examined. 


1. Introduction 


Impulsive and indicial admittances are of great utility in calculating the 
responses of systems to arbitrary applied forces.’ * ° An impulsive admit- 
tance z(t) is the response at time ¢ to the application of unit impulse at time r=0, 
the initial displacement and velocity being zero. For a system having a single degree 
of freedom and subject to an applied force F (t) the displacement is then given by 


t 


0 


where x and dx/dt are zero when t=0. For a system with nm degrees of freedom 
there are n* impulsive admittances 2z,, (t), where this is the response in the r'" degree 
of freedom to the application of unit impulse in the s‘* degree of freedom and the 
initial conditions are that the displacements and velocities are zero. Then, if F, (¢) 
is the applied (generalised) force in the s" degree of freedom, we have 


t 
| (7) Zrs (t 7) dr, . . . (2) 

1d 
where x, and dx,/dt are zero when t=0. The indicial admittances are the responses 
to forces whose graphs are “ unit step functions ” and they are the integrals of the 
corresponding impulsive admittances. Thus the general indicial admittance £,, (#) 

is given by 
t 


0 
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while the general equation corresponding to (2) and obtainable from it by integra- 
tion by parts is 


d 


The expressions given are valid for systems whose dynamical equations are 
linear with constant coefficients. However, an extension to systems which are still 
linear and have time-dependent coefficients can be made ®, but this is not required 
in the present application. 


For the dynamical chains considered in this paper the expressions for the 
admittances are very simple. When the system is “ anchored ” (see Section 2), the 
typical impulsive admittance is 


n 
p=1 


where «, is the p' root of the characteristic equation. The corresponding indicial 
admittance is 


n 
p=1 Wp 
When the system is “ unanchored,” 
Zrs A,t Ais sin . . . . (8) 


p=1 


where w, is the p root of the reduced characteristic equation, and the correspond- 
ing indicial admittance is 


'A 


Ba 22" (1-cos o,f). (9) 
Wp 


o=1 


The coefficient A, is equal to the reciprocal of the sum of the inertias of the 
elements of the system. 


For the kinds of system considered here, and in all cases where the matrices 
of the dynamical coefficients are symmetric, the admittances have the symmetric 
property 


Ars (t)= (t) 
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It may be remarked that the expressions for symmetrically related admittances 
obtained by a given method of reduction commonly appear to differ, but their 
identity can always be proved by making use of the relations between the roots 
and coefficients of the characteristic equation. 


For any given system it will often be best, from the aspect of labour saving, to 
obtain the coefficients A,,, by numerical solution of the basic equations given in 
Section 4. Nevertheless, the explicit expressions for these coefficients given in 
Section 5 are of considerable interest and value. 


A list of integrals which are of use in the calculation of responses by the 
application of relation (2) will be found in Chapter 15 of Ref. 3. The frequency 
equations and receptances of dynamical chains are listed in Ref. 7, while the 
factorisation of the characteristic equations of dynamical chains is considered in 
Ref. 9. 


2. The Systems Considered 


The systems are “ conservative ” and their dynamical equations are linear with 
constant coefficients. In addition, the system is built from elements each of which 
has one degree of freedom, while each element is coupled only to its neighbours on 
both sides and the coupling is purely “ elastic.” Thus we have a chain of elements 
with the adjacent elements coupled through the agency of forces proportional to 
their relative displacement. Some special examples are: — 


(i) Linear system of masses connected by ideally massless springs, the 
masses lying on a straight line and the motion being confined to this line 
(longitudinal motion of mass-spring chain). 


(ii) System of coaxial rigid flywheels connected by ideally massless elastic 
shafts (torsional motion of mass-spring chain). This system is also a 
working approximation to a multiply-cranked piston engine. 


(iii) System of masses connected by ties in tension, the masses lying in one 
straight line when undisturbed, with small parallel displacements of the 
masses in a direction perpendicular to this line (lateral motion of mass- 
spring chain). The tie may be an ideally massless cord in tension, but 
the tension need not be the same in all segments—an example is provided 
by a multiple pendulum hanging vertically under the influence of gravity 
where the tension in any segment is the total weight of the bobs lying 
below it. 


The system is said to be anchored when it is connected to “ earth” through an 
“elastic ” element at one or both ends and it is unanchored when there is no such 
connection to earth. An unanchored system will move uniformly without relative 
displacement of its parts when each of its masses has simultaneously received an 
impulse in proportion to the mass. 
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Whatever be the exact nature of the system, we shall use the displacement 
co-ordinates x,, X.,. . . . Xn, While the corresponding inertias are m,, m.,.. . Mn. 
The “ stiffness ” of the connection between the masses m, and m,4, iS C;,-+,, but the 
stiffness of the connection between m, and earth is denoted by c,,, while c,, is the 
stiffness of the connection between m, and earth. In accordance with the foregoing 
the dynamical equations are 


where F, (t) is the “force” applied to m, at time ¢t. For any free mode of the 
system x, is proportional to exp (i,t), and we shall have accordingly 


X, + C12 — — =0 


— Wp MnJ=0. (14) 


The eliminant of these equations is the characteristic equation of the system. 
This is of the n'" degree in ,* and it is known that the roots are all real and positive 
or zero. However, for an unanchored system the constant term in the characteristic 
equation is zero and, on division by ,’, this leaves an equation of degree ( n— 1) in 
@,*. This will be called the reduced characteristic equation of the unanchored 
system. For an anchored system the roots will, for convenience, be arranged in 
order of magnitude, so that 


As pointed out in Section 3, repeated roots cannot occur. Similarly, for an 


unanchored system 


Corresponding to any non-vanishing root of the characteristic equation there will be 
a modal displacement of the form 


X,=Gry SiN (yf + : se AGED) 
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For the zero roots of the characteristic equation the modal displacement is of the 


form 
X,=at+bt, . ‘ - Gs 


where the coefficients are the same for all values of r. This mode represents a 
uniform motion of the system as a rigid body. 


The expressions for the admittances are valid whichever sign is given to , 
but, by convention and for convenience, all the » will be taken positive. 


3. Some Characteristics of the Systems 


3.1. REPEATED ROOTS OF THE CHARACTERISTIC EQUATION* CANNOT OCCUR 


Since the system is conservative, the occurrence of an r-fold root is associated 
with a mode of motion with (r- 1) degrees of indeterminacy. But in the present 
case there can be no indeterminacy. This follows from equations (12) -(14); for 
let x, and w, be given, then x,/x, is given uniquely by (12). The next equation of 
the set gives 


is + (>) (23 + Cog 
x, Cos x, Coz 


which is unique. Proceeding step by step, it is found that the ratios of all the 
displacements to x, are unique and not infinite. 


Since repeated roots are absent, the discriminant of the characteristic equation 
cannot vanish. This is of some importance as the admittances contain the 
reciprocal of the square root of the discriminant as a factor. 


3.2. THE First AND Last ELEMENTS CANNOT BE AT REST IN ANY FREE 
OSCILLATORY MODE 


If x, is zero, all the displacements must be zero, from the argument given in 
Section 3.1. Similarly, all the displacements must be zero if x, is zero. It can be 
shown likewise that two adjacent masses cannot be at rest in any free oscillatory 
mode. 


3.3. CONDITION THAT A MASS SHALL BE aT A NODE 


Let m, be at a node; it may then be supposed fixed to earth. There are now 
two partial systems, namely, the part of the system on the right of m, and the part 
on its left. Then the characteristic equations of the complete system and of either 
partial system must have a common root and the three eliminants must vanish, but 


*For unanchored systems, the reduced characteristic equation. 


tAt each step the division is by a stiffness coefficient which cannot be zero. Otherwise the 
system would consist of two or more independent portions, 
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this is effectively a single condition. Evidently m, may be given any value without 
affecting the motion in this partcular mode. Hence 


Om, Ju=u, 
where A (w) is the characteristic function and », refers to the mode in question. 


3.4. CONDITION THAT A PAIR OF ADJACENT MASSES SHALL HAVE EQUAL 
DISPLACEMENTS 


In this case the link between the masses is unstrained (for lateral oscillations 
there is a slope node). Let m, and m,., have equal displacements in the mode of 
motion corresponding to ,. Then c,,,, may be given any value without influenc- 
ing the motion in this particular mode. Hence 


Also c,,,4, may be reduced to zero, so the characteristic equation of the complete 
system and of either of the partial systems resulting from cutting the connection 
between m, and m,,, must have a common root. The three eliminants must there- 
fore vanish and this again is effectively a single condition. 


3.5. GENERAL EXPRESSIONS FOR THE PRODUCT AND SUM OF THE CHARACTERISTIC 


Roots 
For an anchored system with n masses the coefficient of w*" in the characteristic 
equation is m,m,.. . .m,, while the constant term is 


where the summation covers the reciprocals of all the stiffnesses. Hence 


Similarly it can be shown that 


CortCie + Cn-1,n+Cno (22) 
m, m, mM, 


Xo" = 


For a singly anchored system (say c,,=0) equation (21) assumes the simple form 


e 
mym,.. 
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ut For an unanchored system with n masses, we have 
Mn 
9) while 
+ . Cos (mz + Ms) + + Mn) 
mm, 
(25) 
\I 
It follows from (21) and (24) that the product of the roots is unaltered when 
. the order of the masses and springs is changed while their values remain fixed. 
of 
‘ 3.6. BOUNDS TO THE VALUES OF THE ROOTS 
It follows from (21) that, for an anchored system, 
c 
Likewise (22) shows that 
te 
m, My, 
Similar inequalities for unanchored systems follow from (24) and (25). 
. Other relations can be obtained by supposing constraints to be introduced, 


which, in general, results in raising the natural frequencies. If, for example, all the 
) connections between the masses are assumed to become rigid, it can be deduced 
- that, for an anchored system, 
Cort Coe) > 

xm (28) 
Exceptionally, this becomes an equality for a symmetric two-mass system. It 
follows from theorems on bordered determinants* that the roots of the character- 
istic equation are separated by the roots of the characteristic equation of the 
system with the first (or the last) mass held fixed. 


) 
3.7. SYMMETRIC SYSTEMS 
The system is symmetric when pairs of masses and pairs of stiffnesses for 
elements symmetrically situated are equal. Thus 
) 
m,=m,, m,=mM,-,. etc. 
(29) 
Co, 
) *See Routh’s Advanced Dynamics of a System of Rigid Bodies, Chapter II. Sixth Edition, 


Macmillan, London, 1905. 
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A mode of motion will be called symmetric when the motions of masses 
correspondingly placed with respect to the two ends of the system are mirror 
images. Thus, for a symmetric mode 


In 
For an antisymmetric mode the relations are 


Suppose first that n is odd. Then the central mass will be at rest in any symmetric 

mode. Hence there are 4(m- 1) symmetric modes. So far as the antisymmetric an 
modes are concerned, the central mass may be supposed divided into equal and re 
separate halves and there are 4(m+1) antisymmetric modes. When n is even, the 
mid-point of the connection between the two central masses will be at rest in any th 
symmetric mode, whereas this connection will be unloaded in any antisymmetric 

mode. Consequently there are 4m symmetric modes and 4n antisymmetric modes. 

All the natural frequencies of the complete system can be obtained from the 
characteristic equations of the simplified symmetric or antisymmetric systems, as 

already indicated. The admittances of symmetric systems satisfy the relations 


Zrs (t)= (1). . (32) 
an 
3.8. SPECIAL SYSTEMS 


For certain special systems all the frequencies and modes can be obtained 
explicitly in simple forms when the number of masses is arbitrary. These systems 
have the characteristics 

and (a) Co) =Cno = 2c, 
or (b) Cor = 2c, Cao = 0, 
Cade, T 


or (c) Co 


Details will be found in Ref. 8. 


4. Calculation of the Admittances 


We take first an anchored system and begin with the calculation of 2,, (f). 
This is the expression for the displacement x, as a function of time ¢ in the free 
motion which follows the application of unit impulse to the mass m, at time r=0. 
The initial conditions for this free motion are accordingly 
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es 
or 
0) 
In accordance with (6) we put 
1) (= Arp Sin 36) 
p=1 
ic 
ic and we suppose that the values of the roots », have already been found. It is 
id required to find the coefficients A, ,». 
le The equations of free motion are (12), (13) and (14). When w, is known, 
y these equations can be solved to yield 
ic 
1S (37) 
x, sp 
) 
and so on. Equations (33) are all satisfied, while (34) and (35) become 
d 


n 
= pA 
p= 


The last equations are linear in the unknowns and can be solved as usual. In this 
way we obtain 2z,,(t). To obtain z,, (t) we use the relations (37), giving 


To find the set of admittances z,,(f), 2..(0), .... 2n.(0), we proceed in a 
manner similar to that already described, but the initial velocities are now 


=... 

(40) 
m, 
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There is in all cases a wide choice of the unknowns. For example, it may be con- 


venient to replace equations (37) by an equivalent set 


(41) 


and so on. Then the unknowns would be the coefficients in the expression for x,. 


For unanchored systems the impulsive admittances are of the form (8) where 


1 


A,= (42) 


That A, has this value follows from the fact that the velocity of the centre of mass 
of the system is equal to the reciprocal of the sum of the masses when unit impulse 
has been applied to any of the masses. The equations corresponding to (38) 
are now 


ay 1 
A,t+ & — 
0 pAiip mM, 
n-1 
A,+ p= 0 


n-1 
A,+ 
p= 


and there are similar equations for the coefficients in the other admittances. 


Labour can be saved by the use of the relations of symmetry (10). Alterna- 
tively, these relations may be used for checking the work. 


5. Explicit Expressions for the Admittances of Systems with One, 
Two or Three Masses 


The impulsive admittances and characteristic equations are here given explicitly 
for both anchored and unanchored systems for values of n ranging from 1 to 3. 
There is no difficulty of principle in proceeding to larger values of n, but the 
expressions become very cumbrous and it is preferable, in any given instance, to 
solve the equations set out in Section 4 numerically. 


Single Mass, Unanchored 


(=, (44) 


1 


10 
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Single Mass, Anchored 
The characteristic equation is 
A (w)=m,0* = (Cy, + =0 
1) 
a, (46) 
Two-Mass System, Unanchored 
” The reduced characteristic equation is 
re 
A (o)=m,m,0? - (mm, +m,)=90 . . (47) 
2) _t sine,t (48) 
my +m, 
8) 
C12 SIN 
Two-Mass System, Anchored 
) A (©) = m,m,o4 [m, Cao) + Me (Coy + + (Cor + C29) + =O. (S11) 


1 
4 = On + Cag — Me) SIN — 
11 (t) M,M,0 ( 12 20 2) 1 
—, Cay — SIN Wt}. (52) 
Cie . 
2 = O= foF Sin ©, sin (53) 


—, (Cho +Coy— SiN wt} . (54) 


Three-Mass System, Unanchored 


The reduced characteristic equation is 


A m,m,m,o* — [c,.m, (m, + m,)+C.,M, (m, + m,)] ©? + 


(mM, +m,+m,)=0.. . (55) 
1 (o,) sin o,f P25 (o,) sin 


where (©) =m,m,o* — + My + Co3)) (57) 
1] 
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An alternative and equivalent expression for the admittance is 


(1) + x 
m,+m,+m,  m,* (m, +m, + - o,”) 


x [—, {((m, + m,) (c,. m,@,") +¢,,.m,} sin o,t+ 


. 


m,o,?=0 
or C,.-—m,,?=0. 


t 
2,,@0= 
uf m,+m,+m, 


Cin" 


+, {(m,+m,) sin of). 


There is a simpler expression for z,, (¢) which, however, fails when 


This failure actually occurs for symmetric systems. The expression is 


The circumstances of failure are discussed in Section 6. 


(t)=2, 


m,+m,+m, mmm, - ©,?) 


Zis ()=2;, (t)= 


m,+m,+m, (," ©,*) 


t 1 


+ 


(t)= 


@,° m,o,”) 


+ 


t [ sin 


[ (Ci. (C23 m,o,”) sin (Ci. M,©,") (C23 — sin Wot 


3 


t Co 


3 3 


12 


x [ (C12 m,,*) sin ot + (c,.—m,,7) SiN Wt 
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where P,, (w)=1,m,04 — [m.c,.+ mM, + ©? + . (66) 

For a symmetric system P,, (o)=P,, (@) 

and we now have also a5, (t)=2,, (ft) 


(t)=2,, (Q=2,, (= 2,5 


Hence there are only four distinct admittances; these may be taken as 2,, (¢), 
(t), 2,5 (t) and 2,, (t). Provided that (c,. m,,*) and (c.; m,©,”) do not vanish, 
we have for any system 


Fg Pg Gr) ‘ . (67) 


This relation may be used to save work in computation or for checking. 


There are alternative expressions for 2,, (tf) corresponding to (58) and (60) 
which are obtainable from these by the substitutions 


m,—> Ms, m,—> m, 


Cig —> C23 —> 


Three-Mass System, Anchored 
A = — (mm, (Cy) + + (Cy. + C23) + MMs (C25 + Of + 


+ + + + MACy, + C12) (Cos +Cg0) + 


+ + +€23Cy1)} — Cy Coy (C4 + C25) — + C39) =0. (68) 


(t) = m,M,M,0,0,0,F [ 0.0, ,”) (w,) sin f+ 
+ 0,0, (w,* ©”) P,; (©,) sin ©,f-+ ©, 0, (w,* ©,”) Sin ] (69) 


and (o)=m,m,o* [m, (C2 + C59) +C23)] (C25 + C59) + (71) 


This is the characteristic function of the system ankylosed by fixing m,. 
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(= 2, (0) 


(,? — (Cos + — M32") SIN Of + ©,7) (Coz + — M35”) SiN ost | 


(H= (W,° SIN ©, f+ 
is ( u(t) M,M,M,,0,0,F L 2) 
+,0, (,?—©,”) sin (@,? — ©,”) sin ost. 
(Q=—— [ 0.0, — + Cig — M,,”) (Coy + — SiN T+ 


+ 0,0, — (Cy, + — 1,02") (Cos + Coo — sin t+ 


+0, — ©,7) (Cy, + — IN, (Caz + — M35”) SiN ost | . (74 


[ 0.0, (@,? — (Cy, + — M1, ,”) Sin 


(t)=2,,(0) = 


+ O30, (o,? m,,") sin of 


+ 0,0. (0,7 (Cy, + sin ost | (75) 


2,, (t)= | 0,0, (w," — P,. (w,) sin + 
33 ( ) M1, 3 3 2 ) is ( if 
+ (W,? P,, (w,) sin + — © ,*) SiN Ost ] (76) 
where 


P,, — [m, + C23) + Me (Coy + C12)] ©? (Cop + (77) 


This is the characteristic function of the system ankylosed by fixing m,. The 
relation (67) again holds good provided that (c,, +¢,.— ,©,”) and (C23 + C39 — m3, 
do not vanish (see Section 6). 


6. Remarks on the Admittances 


The anchored three-mass system will be taken as typical, and the relevant 
formulae are given in equations (68) to (77). The first remark is that F, as defined 
in (70), can be shown to be a form of the square root of the discriminant* of the 


*The reader may be reminded that the discriminant of a function is a form of the eliminant 
of the function and of its first derivative. The vanishing of the discriminant is the condition 
for the occurrence of repeated roots. 
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characteristic equation (68). When the roots are ordered as in (15) F is certainly 
positive (it cannot be zero since repeated roots do not occur, as pointed out in 
Section 3.1). The admittances, which contain the reciprocal of F as a factor, are 
thus never infinite but some of their coefficients may be zero in special 
circumstances. 


It can easily be verified that the general coefficient A,,, contains as a factor 
4,, (©,), where A,,(@) is the co-factor of the (r,s)" element in the characteristic 
determinant. We may go further and give the following explicit general formula 
for the three-mass anchored system : — 


(Wp 427 — p47) Ars 

where ©, is to be interpreted as »,_, when n exceeds 3. It may be noted that 

(p47 — ©4,°) 18 a Square root of the discriminant of the characteristic function of 
the system with the factor (w?— ,”) removed. 


Let w be any root of the characteristic equation A(w)=0. It is found that this 
equation may be written 


P,.. + C39 — - Cas” (Co. + m,o*)=0, 


where P,, (w) is defined by (77). Provided that (c.,+¢;,—m,?) does not vanish, 
this yields 


Cas" (Co; m,w*) 


Cag + Cay — Mw? 


(79) 


provided that c,,+¢C,.—#,? does not vanish. 


always provided that the expressions in the numerator and denominator of (80) are 
not zero. The relation (81) may be used to save computation or as a check on 
the work. 


Equations (79) and (80) are valid for an unanchored system, with c,, and c,, 
made zero, and always provided that the numerator and denominator do not vanish. 
When these equations are used (56) becomes (60). 
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For a symmetric system, 
C23 — 
but for one of the roots the numerator and denominator both vanish and the 
fraction is indeterminate. It can be shown that the true value of the ratio in the 
limit is —1. By use of (79) and (80) simpler forms of (69) and (76) can 
be derived and these are subject to failure in the manner already discussed. 


According to observation, all the coefficients A,,, occurring in the direct 
admittances 2,, (t) are positive, but this has not so far been proved. 
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Compressor-Turbine Matching of 
Two-Spool Turbo-props 


A. W. MORLEY, Ph.D. 


(D. Napier and Son Ltd.) 


Summary: The problem of the best division of pressure ratio between the two 
spools is investigated by means of simplified expressions for component 
characteristics. It is found that a unique position of the L.P. compressor 
operating line (i.e. maximum L.P. flexibility) is obtained when the H.P. 
compressor operating line is vertical. This arrangement requires a low pressure 
ratio on the H.P. spool and gives an engine whose power is sensitive to air 
intake temperature changes, so that the H.P. rotor. tends to overspeed at sea 
level and surge at high altitude. Checks by standard methods show that the 
algebraic approach is useful for understanding the two-spool matching problem. 


1. Introduction 


The “ off-design” performance of a two-spool turbine engine is more difficult 
to predict from the component characteristics than that of the single-spool engine, 
because of the extra labour involved in tracing the exact movements of the operating 
points of the two compressors and two turbines across their respective performance 
maps. To some extent the choice of the leading design data of such an engine 
depends on estimates of how much is lost by “ off-design” operation, so that it 
is important to try to understand the major trends of the matching points at an 
early stage. 


In attempting to forecast “‘ off-design’? performance and as a preliminary to 
more exact work, some use can be made of the algebraic methods described in this 
paper. To be amenable to such treatment the thermodynamic system must be 
simplified. Despite this, checks by more exact methods have proved that deductions 
from the simple algebra are qualitatively correct and provide a better understanding 
of the two-spool turbo-prop performance problem than is obtained by trial and 
error processes among the several characteristics. In particular an interesting point 
is brought to light when considering the best division of pressure ratio between the 
two shafts. 


Received January 1957. 
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L.P H.R H.R L.P 
COMPRESSOR COMPRESSOR TURBINE TURBINE 


FiGurE |. Reference planes. Two-spool turbine engine. 


NOTATION 
A,B,C, constants defined in Section 5 


a=Q/(Ar;) 
k “work-done factor” 
N blade speed (r.p.m.) 
P pressure (total head) 
Q mass flow parameter (0, =W, /7,/P,, and so on) 
r expansion ratio 
L.P. and H.P. compressor ratios respectively (r,=P./P,, ra=P,/P.) 


T temperature (total head) 
t=T,/T, 
W mass flow 


1» %, inlet and outlet air angles 


y ratio of specific heats 


Suffixes 
H,L denote H.P. and L.P. compressor (or turbine) respectively. 
1.2 L.P. or outer compressor 
3 H.P. ori 
HP. engine reference planes (see Fig. 1) 
5.6 L.P. or outer turbine 


Abbreviations 
H.P.,L.P. high pressure, low pressure, respectively 
N.D. non-dimensional. It is, of course, understood that the conventional 
forms of the engine parameters are not strictly non-dimensional 
unless the appropriate linear dimension is introduced, but as this 
paper deals with a fixed configuration the linear dimension may 
be discarded. 
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H.R COMPRESSOR MAP L.R COMPRESSOR MAP 


FiGuRE 2. Compressor maps. 
2. Assumptions 
In the theory the leading assumptions are: — 


Polytropic efficiencies of the compressors are constants over the respective 
performance maps. 


The combustion chamber pressure loss is counted in with the H.P. 
compressor as a loss in H.P. pressure ratio. 


There are no air bleeds, and fuel weight is neglected in comparison with 
the air mass flow. 


The specific heat of the air, and of the turbine gases, is a constant. 


3. Mass Flow Shift between L.P. and H.P. Compressor 


As speed or altitude increases the operating points on the two compressors 
will climb to higher values of N// 7 and the N.D. mass flow corresponding to the 
L.P. operating point will, in the general case, increase. It is a property of the 
two-spool engine that in this type of change the N.D. mass flow increases a greater 
percentage on the L.P. than on the H.P. compressor. Thus, in Fig. 2 the operating 
point on the L.P. moves to the right and farther from the surge line as N,//T, 
increases, while on the H.P. map the movement is nearer the vertical, because the 
change in the N.D. mass flow is less and so the H.P. point moves towards the surge 
line. The reason for this is not far to seek, for it comes from the identity 


If the polytropic efficiency is 7, then r,=P,/P,=(T,/T,)”°~-”, and therefore 
relation (1) may be written in the form 


where §= ‘ ‘ (2) 
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COMPRESSOR EXPONENT 


0-791 


0-7 075. 0-8 0-65. 0:95 
COMPRESSOR EFFICIENCY (POLYTROPIC) 7 


FiGureE 3. Compressor exponent S. 


From relation (la), for corresponding small changes 4Q,, AQ, and Ar, which 
occur with a change in running conditions, it is seen that 


AQ, _ AQ, 3 


As the operating point climbs up the L.P. characteristics with increasing r.p.m. 
or altitude, Q, and r,, increase together. Moreover, S is always positive, depending 
only on compressor efficiency (see Fig. 3). Hence, from (3), the percentage change 
in Q, is greater than the percentage change in Q,, so that the tendency is for the 
L.P. to move out from surge while the H.P. moves towards surge for this type of 
change. 
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On the other hand, if the N.D. L.P. speed is held constant, as when the 
propeller is constant-speeding during a power change in level flight, AQ, and Ar, 
will have opposite signs. From (3), the percentage change in Q, will now be less 
than the percentage change in Q,. In certain cases AQ, may be negative, so that 
the L.P. compressor will move towards surge; in other cases it will be positive, 
and so move away from surge. Under special circumstances which remain to be 
investigated AQ, can be zero, and it is then impossible to surge the L.P. compressor 
by gradual change of power at any altitude. 


4. Changes in the Work-Done Factor 


The “* work-done factor” k, for a compressor, expresses the rise in temperature 
in the compressor in terms of the blade speed or r.p.m., i.e. 


(4) 


To form a mental picture of what changes in work-done factor are likely to 
occur, consider a single-stage wheel with one of the standard forms of blading. 
For example in “‘ symmetrical”? blading the work-done factor is expressed in terms 
of the inlet and outlet air angles z, and z, by the well-known equation 

tan z, — tan 2, 

where C is a numerical constant. If now we imagine the axial velocity decreasing 
while the blade speed is held constant, the inlet angle z, will increase at a greater 
rate than the exit angle z,; in much the same way, when an aerofoil is given an 
increased incidence, the angle of attack increases much faster than the downwash 
angle, and the lift increases until the stall is approached. Hence it is seen from (5) 
that k increases with reduced mass flow at fixed blade speed; i.e. along any 
constant N.D. speed characteristic in the direction of reduced mass flow the work- 
done factor increases as surge is approached. 


In the multi-stage compressor, if the operating point of one stage moves towards 
the surge the others must do likewise. The work-done factor of each stage will 
follow an angle function more or less similar to (5) and therefore k will increase 
as the stage gets nearer the surge. Hence it is concluded that the work-done factors 
of the two compressors will increase if their respective operating points move 
nearer surge, and will decrease if they move away from it. 


With increase of altitude at fixed L.P. and H.P. r.p.m., it has been seen that 
the L.P. moves away from, and the H.P. in towards, surge. The work-done factor 
on the L.P. therefore decreases and on the H.P. increases. As a direct result 
the H.P. turbine inlet temperature must increase, to balance the demand of extra 
work per Ib. mass flow from the H.P. compressor. Meanwhile the work absorbed 
per Ib. of air in the L.P. compressor is decreased, so that a larger part of the L.P. 
turbine output is available for shaft power. 
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5. The Unique Operating Line of the H.P. Compressor 


The free-running inner spool of a two-spool turbo-prop constitutes a 
compressor-turbine equilibrium pair behaving after the manner of the balanced 
pair of a simple turbo-jet. A well-known property of such a pair is that the 
operating point moves along a single line on the compressor map all the while the 
exit N.D. mass flow parameter (Q,) is a constant. This occurs when the flow is 
choking at the downstream end of the H.P. turbine. Elementary algebra will show 
that the operating line is unique and will yield an equation approximating to the 
true equation for this line. 


Consider the matching of the H.P. spool. Let the N.D. mass flow for choking 
conditions at the H.P. turbine inlet to be denoted by A. Then A will be a design 
constant for the engine, i.e. the H.P. turbine design mass flow parameter. Then, under 
all normal running conditions where the inlet to the H.P. turbine is choking, Q,= A. 


The overall rise in the H.P. compressor is assumed to include the combustion 
chamber pressure loss. The identity 


Again, since the inner spool is balanced, the work done in the compressor can 


be equated to the work output of the turbine: — 


(T,-T.) _ 
T;) 


where B is a physical constant. 


Further, if we postulate choking at the L.P. turbine inlet, Q@;=a constant. By 
a similar reasoning to equation (1), Q@, may be written Q,r* or Ar‘, where r is the 
expansion ratio and S$ the simple function of expansion efficiency appropriate to 
the H.P. turbine. Then, since Ar* is constant, r must be constant, and hence the 
temperature ratio across the turbine 7,/7, must be constant. 


=C,aconstant, . . (8) 


Therefore 
where C depends only on the design pressure ratio (or heat drop ratio) of the 
inner spool turbine. 


By combining (6), (7) and (8) and substituting for 7,/7., the equivalent value 
ry°°-*, the following law for the H.P. compressor operating line is obtained: — 


It is readily seen that (9) represents a single curve (for a given compressor 
polytropic efficiency) across the H.P. compressor map. A is a constant fixed by 
engine size, B is principally a physical constant depending on the specific heat of 
the working fluid and C is the chosen value of the H.P. turbine work parameter. 
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Equation (9) holds for change of altitude and change of r.p.m.. As for a 
simple jet engine, if the work-done factor of the compressor is constant, the H.P. 
turbine inlet temperature is constant when the H.P. r.p.m. is constant, a result 
which follows on combining (4), (7) and (8). Thus the law expressed by (9) contains 
the special case of constant turbine inlet temperature and constant H.P. r.p.m. if the 
H.P. compressor work-done factor is constant. For this to be so the operating 
line would require to run roughly parallel to the surge line. The law (9) applies 
only when both H.P. and L.P. turbines are choking at their respective inlet planes. 


6. The L.P. Compressor Operating Point 


It has been seen that any changes in Q. will be magnified on the L.P. compressor 
map if they are accompanied by an increase in L.P. pressure ratio. This occurs 
with increase of altitude or propeller speed. Since it is desirable to resist the 
movement of the L.P. point away from the region of good compressor performance, 
it is necessary to limit the change in Q,. From (9) it is seen that the variation 
in Q, over a given range of ry is proportional to A /(BC). For a fixed size of 
engine, A and B are determined by the design. The remaining constant C can be 
reduced by dropping the design value of ry. But when a definite overall design 
pressure ratio is required a reduction in design ry requires an increase in design r,. 
It is not difficult to see that Q, would increase again if r, is increased, which would 
offset the advantage gained by restricting the change in Q,,. 
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DESIGN TEMPERATURE RISE (°c) 


FicurE 4. Percentage change in pressure ratio for intake temperature change from 288°K to 
216°K. Calculated from equation (10) for constant polytropic efficiency=0°85. 
T, (mean)=247°K., AT, /T,=0-29. 
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Thus we may write down a compressor equation from (2) and (3) in the form 


r2a-s) =l+kin,7, . ‘ (10) 


where n, is the N.D. r.p.m., N,/ /7,; and from this it can be shown, as in Fig. 4, 
that for a typical change in engine intake temperature, AT, /7, =0-29, the percentage 
change in r;, increases almost in proportion to the design compressor temperature 
rise. Therefore the greater the design pressure ratio the greater is the percentage 
change in r, and the greater the movement of the L.P. operating point for a given 
intake temperature change. Qualitatively a similar result applies to the H.P. 
compressor. At constant r.p.m. both suffer similar changes in intake temperature 
over the altitude range. 


To follow the combined effect of pressure ratio changes on the L.P. operating 
point, it is first necessary to introduce a new function to define the fuel charge or 
throttle position. A convenient parameter is the overall temperature ratio between 
H.P. turbine inlet and L.P. compressor inlet. 


Let T/T, =¢. 
The identity 0,/Q,=(7,/7,)“P./P,) (P,/P.) can now be written in the form 


40, _ Orn _ 1 


In other words, the percentage change in L.P. N.D. mass flow is always the 
sum of the percentage changes in the L.P. and H.P. pressure ratios less one half 
the percentage change in overall temperature ratio. 


or, for small changes, 


Two types of change will be examined in the light of equation (12): first, a 
change of altitude at constant r.p.m. and second, a change in power at a fixed 
altitude. 


7. L.P. Compressor. Altitude Change at Constant R.P.M. 


To examine the effect of intake temperature change on the movement of the 
L.P. point, relation (12) is first converted into a more convenient form. With 
constant polytropic efficiency, the overall pressure ratio may be expressed in terms 
of the overall temperature rise (7,/7,) thus: — 


Assuming constant work factors on both compressors, (7;—7,) and 7, 
are constants, in which case the differential A{(T7,—7,)/T,} is simply 
—(T, -T,)(QAT,/T,?) and At/t is - AT,/T,. Substituting into (12) gives the result 


4Q,/Q, = = T)/T, (13) 
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Equation (13) is for a fixed maximum temperature 7, and does not contain the 
heating parameter. 


A somewhat similar equation is obtained for the change in L.P. pressure ratio, 
from (10), namely: — 


AT,/T, 
Curves calculated from (13) and (14) are plotted in Fig. 5. The polytropic 
efficiency is assumed to be 85 per cent. 
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INLET AIR TEMPERATURE (‘k) 


FiGurE 5. Percentage changes of non-dimensional mass flow Q, and L.P. pressure ratio with 
percentage change in inlet temperature. Calculated from equations (13) and (14) for constant 
polytropic efficiency =0°85. 


A AT T 
2, 1 (equation (13)). ————— s—+ / meal (equation (14)). 
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In Fig. 5 it will be noted that the quantity given by (14) has been multiplied 
by S before plotting alongside AQ,7,/(AT,Q,). This is to permit examination in 
the light of equation (3), which shows that, if the change in AQ, were made zero, 
then 


SAr;, _ AQ. 
ry oF 


By interpolating between the three values for L.P. temperature rise it appears 
that A4Q/Q, can be made sensibly the same as SAr,/r;, and so to satisfy the last 
equation, if the L.P. temperature rise is about 210°C on the 14:1 engine and 
about 185°C on the 12: 1 engine. The L.P. pressure ratios would then be 5-1 and 
4-4; the H.P. pressure ratio is found to be 1:73 in both cases. Thus it appears 
possible, within the stated assumptions, to divide the full pressure ratio between 
the two spools so that over a wide range of intake temperature there will be no 
change in N.D. mass flow on the H.P. compressor. 


The condition for zero shift on the H.P. compressor map can be found from 
(13) and (14), equating SAr,/r;, and AQ,/Q, for equal intake temperature changes. 


This gives the result 


: 
S S 
This relationship connecting the overall rise and L.P. rise with intake tempera- 
ture is seen to be equivalent to 


T,/T,=8. 


The solution given by (15) and (16) is illustrated in Fig. 6. The special value 
for the H.P. pressure ratio is quite independent of the overall pressure ratio. 


The same result is found on calculating the stationary (minimum) value of Q. 
from (9). Since A, B and C are constants, Q, is stationary when the rate of change 
of (1/ra** 1/rn*) with respect to ry is zero, i.e. when 


Forms (16) and (17) are seen to be consistent. With 85 per cent polytropic efficiency 
rx is 1-728, with 95 per cent it is 1-719. 


With this special value of H.P. pressure ratio the change in Q, is zero, i.e. the 
H.P. operating line is vertical. As altitude increases and 7, falls, the H.P. turbine 
N.D. speed must increase so that the vertically climbing operating point moves fairly 
rapidly towards surge. The work-done factor then increases and one would expect 
a corresponding increase in 7,, the turbine admission temperature. In arriving 
at (13), and again (17), a constant maximum temperature has been assumed as the 
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FicurRE 6. Compressor pressure ratios for constant non-dimensional mass flow on H.P. 
compressor. Calculated from equations (15) and (16) for constant polytropic efficiency=0°85. 
A for 14:1 overall in stratosphere, B for 12:1 overall in stratosphere. 


desirable condition for the free running H.P. spool at constant rotational speed. 
It appears that this assumption can only be approximately true for the two-spool 
engine and that generally an increase in maximum temperature will occur as altitude 
increases. The engine with constant Q, (i.e. ra= 1:73) appears likely to have more 
top temperature increase with altitude than one with increasing Q, (i.e. rx > 1-73), 
since the former moves more sharply towards the high altitude surge boundary. 


Some further examination of this point is made later, in Section 9. 


Figure 5 is of some interest for comparing changes in L.P. mass flow. The 
operating point moves away from surge on the L.P. compressor as the intake 
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temperature falls. It is seen that the mass flow shift is increased (numerically) if 
the overall pressure ratio is increased. Thus, changing from 12:1 to 14:1 
overall gives the same effect as an extra 20 — 30°C of inlet temperature over which 
the L.P. compressor must operate satisfactorily. This does not depend upon the 
division of pressure ratio between the spools. 


8. L.P. Compressor Operating Line. Effect of Power Change at 
Constant Altitude 


We turn now to the case of varying fuel flow at a fixed altitude and constant 
propeller r.p.m. The L.P. point will be located on a line of constant N,//7,. 
If the running point moves only slightly along this line over the throttle range, so 
much the better from the point of view of L.P. flexibility. 


With Q, and r,, constant, equation (12) reduces to 


Thus for a fixed L.P. operating point we should require 


where a is a constant (=Q,/(Ar,)) for the particular altitude and L.P. speed. In 
other words, over the power range at constant L.P. speed the heating parameter 
T,/T, must vary as the square of the H.P. pressure ratio, for a fixed L.P. point. 


The equation of the H.P. operating line (9) may be transformed on to the L.P. 
operating map by making use of (1) and (11), which hold irrespective of the operat- 
ing conditions. If we may assume that both compressors have the same polytropic 
efficiency we obtain, for the L.P. operating carpet, 


Figure 7 shows some results calculated from (20) for a polytropic efficiency of 
85 per cent (S=0-832), for an engine with a design overall pressure ratio in the 
stratosphere of 14:1. The lines shown give the pressure ratio and mass flow 
parameter for the L.P. compressor for two values, 4 and 4:5, of the heating 
parameter ¢. (At a given altitude higher values of t give increased power output.) 
In the upper diagram the H.P. compressor has a design rise of 104°C and it is 
seen that an increased power brings the operating line nearer the surge. In the 
lower diagram the H.P. compressor has a nominal rise of 208°C and here with 
increased power the movement is away from surge. 


Consider a power change at fixed altitude in the light of equation (20). Since 
A, B and C are design constants, the net effect on the N.D. mass flow Q, of a 
change in heating ratio t depends only on the value of the L.P. pressure ratio r, 
at the instant. 
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If the shift in Q, is to be small, in which case the change in r, will also be small, 
(20) may be differentiated with respect to ¢, treating Q, and r; as constants. By 
making use of (11) we obtain 


as the condition which would allow the L.P. operating point to remain fixed. Now, 
from the definitions of B, C and rf, 


whence, from (21) and (22), 


Thus we arrive again at the special value of the H.P. pressure ratio (found in 
(16) and (17) when considering intake temperature change), giving zero mass flow 
shift on the H.P. compressor characteristic, i.e. a vertical equilibrium running line 
for the inner spool. 


The simple algebraic solution given suggests that this value of rx for maximum 
flexibility depends solely on compressor efficiency. In an actual engine ry will 
increase with H.P. r.p.m. and it will only be possible to approach a constant value 
if the H.P. r.p.m. band is relatively narrow. Suppose ry is chosen to suit the 
average cruising H.P. r.p.m.; then at the ends of the operating range the L.P. 
running point will move relative to the surge line as shown in Table I. 


TABLE I 
Rati Ch Movement of L.P. 
H.P. r.p.m. Operating Point at 
Constant L.P.r.p.m. 
{ Increase Towards surge 
Below mean <1°73 | Decrease Away from surge 
{ Increase {No movement 
Mean 173 | Decrease | (unique line) 
ar { Increase Away from surge 
Above mean tiles | Decrease Towards surge 


The critical condition of approach to surge for the L.P. compressor is low 
N.D. speed—low altitude, and for the H.P. compressor, high N.D. speed—high 
altitude. The upward movement of the H.P. point towards surge at altitude causes 
an increased turbine inlet temperature, i.e. the engine must be temperature limited 
at altitude. 


At sea level take-off conditions, on hot days, with the L.P. point nearest to 
surge and the H.P. point farthest from surge, the turbine inlet temperature will fall 
well below design and power will be lost. If more fuel is supplied to restore the 
turbine inlet temperature, the H.P. spool will overspeed by an amount which will 
be greater the hotter the day. If the L.P. (i.e. propeller) shaft is allowed to over- 
speed, T, increases, and therefore Nu//T, decreases, so that the H.P. point 
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FicurE 8. Data on H.P. and L.P. turbine maps (diagrammatic). Two-spool propeller turbine. 
The line AB represents movement of operating point from sea level static to altitude maximum 
cruise conditions. 


moves still farther from surge, which again lowers turbine inlet temperature and 
decreases the power. Conversely, if L.P. speed is reduced it will be possible to 
increase the fuel without exceeding limiting H.P. speed, and this may be sufficient 
to offset the other effects of reduced L.P. speed and give increased power. It would 
therefore be expected that the power curve at a fixed altitude would rise sharply 
with H.P. speed but would be fairly flat over a range of L.P. speed. 


9. The H.P. Turbine 


The turbine performance map gives the data shown in Fig. 8. It consists of 
a series of roughly parallel constant N.D. speed lines, falling nearly vertically as 
the load parameter (Tn -— Tou:)/Tin decreases from the limiting loading line towards 
the horizontal WN /P axis. 


It has been seen that the H.P. turbine in the normal operating range runs at a 
constant temperature ratio 


T, 
If the H.P. operating point moves at all it must move parallel to the WN /P axis. 


The magnitude of the movement: depends on the change in work-done factor ky 
of the H.P. compressor according to the equations (4), (7) and (8), i.e. 


7,-T. BC 
23 
5) 

As already shown, ky is greatest at altitude and least at maximum L.P. r.p.m. 
at sea level. Given a constant H.P. r.p.m., equation (23) shows that 7, will increase 
or decrease proportionately to ky, and if the H.P. r.p.m. increases 7, will increase 
as the square. 


The effects just noted are responsible for the sensitivity of the turbine inlet 
temperature to changes in engine intake temperature. In a design study using actual 
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compressor characteristics, ky fell by 6 per cent between 36,000 ft. and sea level, 
while the turbine inlet temperature dropped 58°C (5 per cent of its absolute value). 
On a hot day (1.S.A.C. +30°C) at sea level, ky fell a further 34 per cent and the 
turbine inlet temperature dropped another 40°C or 34 per cent. In a second 
example, with the pressure ratio divided so as to give a vertical H.P. characteristic, 
ky fell 11 per cent between 36,000 ft. and sea level, with a turbine temperature 
drop of 138°C or 12 per cent. The power output of the engine will, of course, 
follow the turbine temperature and the losses in potential output will be magnified 
as ky moves from its area of maximum values. The degree of the temperature 
sensitivity is thus greater for the engine with the steep H.P. operating line, which 
must react against the good flexibility of such a design. 


The H.P. outlet temperature is given by 


from which it is seen that 7, varies only with ky at constant H.P. speed. Then, 
since ky decreases, both the turbine inlet temperatures decrease when the engine 
inlet temperature rises. 


It is seen that the off-design movement of the H.P. turbine point will be small. 
This will simplify the aerodynamic design of the turbine and ensure that its efficiency 
does not change appreciably over the normal running range. 


10. The L.P. Turbine 


Since the L.P. turbine exit does not remain choked, there is no condition 
corresponding to (8) to govern the work parameter of this turbine. There will also 
be a change in WN,/P, to consider. 


With increase of altitude, at constant H.P. and L.P. speeds, it has been seen 
that the increase in the work-done factor ky will increase T, (see equation (24)), 
and will therefore reduce the L.P. turbine speed parameter N,,/ /7,. The percentage 
decrease in N,/ /7, will be half the increase in ky. The L.P. turbine point thus 
moves left with respect to the speed lines. As the engine climbs the L.P. turbine 
operates nearer the choking condition, so that the point also moves upwards with 
respect to the load lines. With increased fuel at a fixed altitude and L.P. r.p.m., 
T, increases along with Ny and 7, (see equations (23) and (24)) and so the N.D. 
speed of the L.P. turbine decreases as power increases. 


The limiting load line of a turbine falls at the lower N.D. speeds. Thus the 
L.P. turbine most nearly approaches the limiting load condition at maximum 
altitude, maximum power. 
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COMPRESSOR-TURBINE MATCHING 


A decrease in the work-done factor k, of the L.P. compressor occurs at 
altitude and is equivalent to a reduction in L.P. compressor work which appears 
as a gain in L.P. turbine output. 


In all the foregoing it has been assumed that the propeller power is taken off 
the L.P. shaft. This permits balanced running of the inner spool, upon which most 
of the properties discussed will depend. 


11. Conclusions 


Although the conclusions to be drawn are only of a qualitative nature, this 
application of simple algebra to the two-spool matching problem helps to clarify 
some of the complexity of the engine performance and assists in the choice of the 
division of pressure ratio. 


By choosing the H.P. pressure ratio at about 1-73, independent of the overall 
pressure ratio, a unique operating line is obtained on the L.P. compressor map 
as well as on the H.P. compressor map, irrespective of the turbine inlet temperature. 
(Note: The H.P. pressure ratio includes the combustion chamber pressure loss.) 


When the H.P. pressure ratio exceeds about 1-73 the L.P. compressor tends 
to run into surge if the power is reduced at constant L.P. speed, and when it is 
less than 1-73 it runs away from surge. 


If the H.P. pressure ratio is about 1-73, the H.P. operating line is vertical, 
going into surge at high power, high altitude. If it exceeds 1-73, the H.P. operating 
line moves towards surge but less rapidly. Due to the approach towards surge the 
turbine inlet temperature increases with increase in N.D. H.P. speed. The two- 
spool turbo-prop is, in consequence, r.p.m.-limited at sea level, i.e. temperatures 
must be lower there than the design maxima, in order to avoid overspeeding. At 
high altitude the engine is temperature-limited, i.e. H.P. speed must fall below the 
maximum if the design turbine inlet temperature is not to be exceeded. 


The two-spool turbo-prop, at constant r.p.m., is sensitive to intake temperature 
change. The most flexible arrangement (rg=1:73 approximately) is the most 
sensitive. 
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The Application of a Block Diagram 
Representation to Aircraft Stability Analysis 


F. A. SUMMERLIN and N. SULLIVAN 


(Sperry Gyroscope Company) 


Summary: The block diagram representing an aircraft in either its lateral or 
longitudinal motion is given and is used to explain the aircraft response in the 
antisymmetric stick-fixed case. A generalised autopilot is added to the aircraft 
diagram and the method of analysing the combination by conventional servo 
techniques is indicated. It is shown how an estimate can be made of the effects 
on stability of non-linearities in the output member of the control system. 


1. Introduction 


The object of this paper is to develop a method of representing an aircraft, or 
an aircraft-autopilot combination, in such a manner that techniques for the analysis 
of feedback systems, which have become highly developed in the past decade, can 


be readily applied. 


Before attempting to represent the aircraft in this form it is of interest to 
indicate the general form of the equations of motion of the aircraft and the 
conventional methods of dealing with them"). 


NOTATION 

X,Y,Z 

a, B,y 

F (8), F (5), F3(s) 

6 

Ww 

u 

6, 


Laplace operator 

general aerodynamic variables 

general complex aerodynamic derivatives 
general aerodynamic forcing functions 
incremental pitch angle 

incremental normal velocity 

incremental forward velocity 


general input function 


Received January 1957. 
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6, general output function 
K general gain function, real and independent of frequency 
G(s) general gain function, frequency variant, complex 
M sgeneral aerodynamic moment (or force) 
@ incremental roll angle 
w~ incremental yaw angle 
v_ sideslip velocity 
relative density 
i, moment of inertia in roll 
i, moment of inertia in yaw 
i, product of inertia 
Mr, moment/angular velocity derivatives 
l,,n, moment/velocity derivatives 
lift coefficient (=4C,) 
k’=k, tan 0 
yp» force/angular velocity derivatives 
y, force/velocity derivative 
R,,R,,R, stability roots of the simple roll loop 
Y,,Y.,Y, stability roots of the simple yaw loop 
X,,X,,X, zeros of the roll/ yaw forward feed 
Z,,Z, zeros of the yaw/roll forward feed 
C,,C,,C;,C,,C; stability roots of the outer loop 
A,,A,; short term automatic control functions 
O,,O, power output member characteristics 
A,,A, long term automatic control functions 
A,,A, auxiliary automatic control functions 


2. The Equations of Motion of an Aircraft 


The stability of an aircraft, for small deviations from steady rectilinear flight, 
can be calculated from a knowledge of the aerodynamic derivatives by an application 
of the elementary principles of the dynamics of rigid bodies, employing moving axes 
of reference defined by Euler angles. 


Briefly, the aerodynamic, gravitational and inertial forces and moments are 
equated to zero, giving six simltaneous equations, which, in the absence of 
substantial cross-coupling terms, divide into two separate groups of three linear 
differential equations each. 
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The methods of the operational calculus can be used to write either of these 
groups in the form 


and to solve them algebraically. In this array, x, y, z are the physical variables, 
and the coefficients 2, 8, y, etc., are functions of s, the complex Laplace 
operator. F,(s), etc., are the Laplace transforms of applied forces and 
moments. 


The full solution of equations (1), (2) and (3) is 


F,(s) [2 F (s) a, B, F,(s) a, B, 
F,(s) B; Ys a, F;,(s) a, F,(s)| 


The response in x, to a forcing function F,(s) only, is therefore 


B, Y2 
|%3 Bs Ys 


and hence the transfer function of the system in this mode can be defined thus 


B, 
By 


The corresponding transfer functions for other modes can be similarly written 
down by inspection. 


While this procedure is both facile and exact, it gives little insight into the 
interaction of the various aerodynamic derivatives. The dominant relationships 
necessary for satisfactory stick-fixed stability can be traced from the factors of the 
characteristic equation’ through its coefficients, back to the terms of the stability 
determinant which are identified with the aerodynamic derivatives (Ref. 1, pp. 
94-177). This provides simple approximate rules for the design of acceptably stable 
aircraft, but to pursue this technique into the study of automatic control and other 
problems of increased complexity involves considerable labour. 
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BLOCK DIAGRAM REPRESENTATION 


38. The Aircraft as a Feedback System 


The block diagram is a basic tool of the engineer engaged in the analysis of 
feedback systems. It is a plan of the system which shows the relations and 
interconnections between the various physical variables. Each functional component 
of the system has a “transfer function’ which is the Laplace transform of the 
differential equation relating the output of the component to its input. It is 
represented on the plan as a block in which its transfer function is written. The 
various blocks are joined by lines bearing arrows indicating the direction of flow of 
information. A number of flow lines can only converge when they represent the 
same sort of thing, and then the outputs of their originating blocks are to be added 
together before being multiplied by the transfer function of the subsequent block(s). 
Where a flow line divides, each branch carries the whole output of the 
previous block. 


Instead of the more or less complex array of components of most servo- 
mechanisms, each with a relatively simple transfer function, the aircraft is a single 
entity with a complicated system of transfer functions. The block diagram is 
therefore most readily derived from the equations of motion, (1), (2) and (3), the 
symmetric and antisymmetric modes still being treated independently. Each of the 
three variables, #, w and u in the symmetric case, and 9, ¥ and v in the antisymmetric 
case, is isolated from the equation that it dominates, e.g. ¢ is obtained in terms of 
the other variables from the equation of rolling moments, giving, in the general case, 


Each of these equations can be represented by a block diagram, and these 
are shown in Fig. 1. The equations are simultaneous, so that the blocks must be 
joined together: as‘ in Fig. 2, which is a complete representation of the dynamics of 
the symmetric or antisymmetric motion of the aircraft. 


Figure 2 is the diagram that has to be drawn before the equations of motion 
can be set up on an analogue computer. Although a computer enables the aircraft 
stability to be studied in an empirical manner, and the effects of the individual 
derivatives can be estimated by changing their values one by one, it gives no 
insight into the underlying theory. Although this diagram demonstrates clearly 
the symmetry of the equations from which it is derived, it needs to be redrawn to 
give Fig. 3, from which it is easy to derive a form more suitable for feedback 
analysis. This can only be done, however, by splitting one of the blocks 1/2,, 1/8, 
or 1/y,, with the result that one of the variables and the corresponding forcing 
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Ficure 1. The block diagrams of the three equations of motion of one aircraft stability mode. 


z 
~B, 
“Bs y 


FicurE 2. The block diagram of one aircraft stability mode. 
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Rearranged block diagram of one aircraft stability mode. 


Ficure 4. The simplified block diagram of one aircraft stability mode. 
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(8) 


FEEDBACK 
BLOCK 


FORWARD 
BLOCK 
FiGurE 5. The general feedback system block diagram. 


function do not appear explicitly. In this case the diagram has been redrawn so 
as to eliminate variable z by splitting the block 1/y,. Thus we have four 
blocks to add to the remaining parts of the diagram, to give the less complicated 
diagram shown in Fig. 4. 


These are: — 


(a) Into the block 1/2, 


[ + 2,1 |x, a feedback across 1/2, which forms a loop which will be 
Ys called loop 1, 


and [ + Bs uy y, which joins in parallel with the cross-coupling block - 8, 
Ys to form a forward-feed block, known as feed-forward 4. 


(b) Into the block 1/8, 


[ + Bs "2 y, a feedback across 1/8, forming loop 3 


3 


and [ + a, 22] x, which, with —z,, forms the feed-forward block 2. 


The input and output of 2 feedback loop are related by the well known 
expression 
_ K,G,(s) 
6, 1 K,G,(s)K.,G.(s) 


(9) 


where 4, is the input function, 6, is the output function, and K,G,(s) and K,G.(s) 
are the transfer functions of the forward and feedback parts of the loop respectively 
(Fig. 5). The product [K,G,(s)] [K.G.(s)] is called the “‘open loop transfer function” 
(Ref. 4, Chap. 4). The negative sign in the denominator of expression (9) arises 
because it has been assumed that an addition takes place at the input of block 
K,G,(s) (i.e. positive feedback). 
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| * % B, 


B3y,-B, ¥3 
% 


Ficure 6. Further simplification of the block diagram of one aircraft stability mode. 


This relationship is here used, first to further simplify the two loops, giving, 
as shown in Fig. 6, 


1/2, Ys 


for loop | 


Ys 


“a 


Boys 


for loop 3, 


and it is then applied to the system to determine the desired transfer function. 
For example, 


F (s) 7 (2372 - 
— Bs¥2) (2:¥s 


¥3 (B27; = 
— Bsy2) — 371) (2372 — 22Ys) 


(Bays Bs 2) (10) 
[41 (Bovs — Bs¥2) + By +71 (4283 — 


When the factor y, is cancelled, equation (10) is the expansion of equation (5), 
thus confirming the analysis. The fact that the denominator of whichever block is 
split cancels from transfer functions such as x/{F,(s)}, x/{F.(s)}. y/{F,(s)} and 
y/{F.(s)} is due to the symmetry of the block diagram, which arises from the 
symmetry of the equations of motion. 
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The block diagram can be used, however, to determine other quantities of 
interest besides the variables x, y and z. For example, it is a simple matter to 
determine the transfer function M(y)/{F,(s)} where M(y) is the aerodynamic 
disturbance at the input to the block 1/8, due to a disturbance F,(s) applied to the 
input to the block 1/2,. (In an actual case F,(s) could be an applied disturbance 
in roll and M(y) the resultant aerodynamic yawing moment.) In transfer functions 
of this kind the factor y, remains in the denominator as one of the stability roots. 


4. The Application of the Diagram to Antisymmetric Stick-Fixed 
Stability 


4.1. THE EXAMPLE CONSIDERED 


The block diagram of Fig. 4 is a pictorial representation of the interactions of 
moments and forces which are responsible for the inherent stability characteristics 
of the complete aircraft. The component parts of the diagram are comparatively 
simple and if, for instance, we consider the stability of the aircraft with controls 
fixed, it is possible to analyse the aircraft performance in easy stages, so that 
particular characteristics of the complete aircraft can be traced back to particular 
characteristics of the individual blocks, and hence to the derivatives of the aircraft. 
This is an important step in that it makes the correction of any undesirable response, 
either by changes in aerodynamics or by the introduction of automatic stabilisers, 
a much more logical process. 


As an example, this section is concerned with the antisymmetric stability of 
an aircraft and it is shown how the Dutch roll oscillation originates in one particular 
loop, and how control terms suitable for automatic stabilisation can be derived, 
virtually by inspection of that part of the diagram. The actual diagrams correspond 
to a particular aircraft at low speed and altitude. 


In the equations of motion for this case the three variables are roll angle 9, 
yaw angle y, and sideslip velocity v. Since the control moments applied by 
movement of the ailerons and the rudder affect principally ¢ and ¥ respectively, 
it is generally convenient to split the block that develops the remaining variable, v. 
Using the non-dimensional form of the derivatives'?) the block diagram becomes 
as Fig. 7. 


4.2. Loop 1: StmpLeE Loop 


This loop describes the aircraft motion in pure roll; i.e. it is assumed that no 
yawing motion occurs. The forward part of this loop, p./{s(si,—1,)}, is derived 
from the moment of inertia and damping of the aircraft about its roll axis. The 
feedback part of the loop, /, (sy, + 2k1)/ {ms (s— yy)}, represents the rolling moment 
due to the sideslip velocity produced by the displacement of the aircraft in roll. 
The open loop transfer function is thus 


S+ (2 /Yp) ky 
l E [s — (1, /i.)] [s - 
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s(si, — Hy s(si.- 
2 
#, (8 -y v) 


4 
(si,+4,.) _46 [He K’) 
He 


FicurE 7. The block diagram of the antisymmetric stability mode. 


and the closed loop transfer function becomes, by equation (9), 
{s(sia 


Pe (s aa Vv) (1 1) 


which reduces to 


Figure 8 illustrates how the roots of the equation 


is [s (1, /is)] [s 


[i.e. the poles of the closed loop transfer function (11)] move, as the part of the 
open loop transfer function independent of s is increased from zero to its actual 
value y,l,/i,. It therefore indicates how the poles of the closed loop transfer function 
are produced from the singularities of the open loop transfer function, under the 
influence of the feedback.* 


*Loci of this type were first suggested by Evans‘®) and are a powerful graphical method of 
obtaining the stability roots of a system. Simple rules have been devised for the rapid 
sketching of loci, one rule being that the loci start at the poles of the open loop transfer 
function at zero loop gain and terminate at the zeros at infinite loop gain. The form of the 
loci depend only on the singularities of the open loop transfer function and the position of 
the closed loop poles on the loci only on the loop gain. Refs. 4 and 13 also give 
descriptions of the technique. 
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OPEN LOOP TRANSFER FUNCTION 


y 


ZERO AT~38-8 NOT SHOWN ON DIAGRAM 


\ 


FicureE 8. Stability root diagram, simple roll loop. 


The closed loop poles include a rapid subsidence R,, which originates from, 
and is not far removed from the point /,/i,, i.e. the influence of the sideslip feedback 
on this mode of motion is relatively slight, and the motion is approximated to by 
considering only the inertia and damping in roll. The substantially oscillatory pair 
of poles R, and R,, originating from the open loop poles at y, and the origin, are 
generally not of importance in the stick-fixed case, since it will be seen on closing 
the outer loop that they return to the real axis. They do, however, represent a 
mode of motion which would occur in an aircraft tightly controlled in yaw by an 
autopilot, but uncontrolled in roll. 


Having determined the roots of the denominator, expression (11) may now be 
rewritten. 


(uz /is) (5 Yo) 
(s— R,)(s— R,) (s- R;) 


(11a) 


4.3. Loop 3: SimpLE Yaw Loop 


This loop describes the motion of the aircraft in yaw, assuming that no rolling 
motion occurs. The forward part of the loop, »./{s(si--n,)}, is derived from the 
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inertia and damping in yaw, while the feedback part of the loop, 
— S Yr) Hok’}/ [He (s— 


represents the yawing moment due to the sideslip velocity produced by the displace- 
ment of the aircraft in yaw. 


The open loop transfer function is therefore 


S (si, N,) (S — yx) 


and the closed loop transfer function is, by equation (9), 


(S— yx) 
— + Nr) +S + My (Me — Yr)] + 


(12) 


In general the closed loop poles comprise the complex pair, Y, and Y,, and 
Y,, which is on the real axis close to the origin. In level flight k’ is zero and Y, 
is at the origin, expressing the indeterminacy of the aircraft heading. In climbing 
or descending flight, the component of gravity along the X-axis of the aircraft moves 
this pole away from the origin, rendering this neutral stability divergent in climbing 
flight and convergent in descending flight. 


The principal simple yaw poles, Y, and Y,, can be evaluated explicitly, in level 
flight when k’ is zero, as 


le l r 3° JF 


which, since n,, y, and y, are usually small compared to »,, reduces approximately 
to Ly, +(n,/i.)] /ic). 


Expression (12) can be rewritten for the level flight case as 


(42 /ic) (S Ww) 


4.4. FORWARD-FEED 2: ROLL TO YAW 


This block introduces the direct interaction terms i, and n,, which describe 
the yawing moments caused by rolling accelerations and velocities, together with 
the interaction via slideslip velocity. The transfer function is 


Ny (SYp + 
Ho (S— yr) 


which reduces to 


Sie + (n, leYr) +S (Np Nyyv) bok Ny 
(S— 
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The numerator can be factorised to give, say, 


i, (s— X,) (s— X,) (s X;) 


that is, the expression has a pole at y, and three zeros X,, X, and X;. 


4.5. FFORWARD-FEED 4: YAW TO ROLL 


This block, describing the rolling moments due to yawing accelerations and 
velocities, both direct and via sideslip velocity, has the transfer function 


Sig +sly [5{1— + KI 
which reduces to 
Si, +8? + Poly) Mok’ly (15) 
In level flight k’ is zero, and hence equation (15) becomes 
which can be factorised to give rm | 
e 
; elim 
My (Ss — yy) 
of fr 
4.6. THE COMPLETE AIRCRAFT Loop upo. 
The open loop transfer function of the outer loop of Fig. 7, which represents ae 
completely the dynamics of the aircraft for small perturbations from rectilinear level oe 
flight can now be obtained by multiplying together expressions (lla), (12a), (14a) thta 
and (15a) giving the 
(16) 
lal (s—R,) (s— R,) (s— R;) (s- Y,) (s—Y,) app 
This expression can be equated to unity to give a characteristic equation of 
the antisymmetric stability of the aircraft. The poles and zeros of the function (16) 
are shown in Fig. 9, together with the stability roots, C,,. . . ., C;, obtained by the 
methods given in Refs. 6, 4 and 13. whi 
The root C, has the value of y,. This must occur since reference to equation 
(10) and later paragraphs shows that a factor y, (which is equivalent to s— y, in this s(s. 


block diagram) always cancels with a corresponding zero. The remaining roots, 
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OPEN LOOP TRANSFER sf C 
FUNCTION L 
0000065 Y; 
(s+ 3-9)(s*+0-18+1-7) (s*+0:28+6-9) 
4” 3 2 2 3 4 
ZEROS AT ZERO aT 
R, (s) 
2- 
3 NS; 
FIGURE 9. The complete stability root diagram of the antisymmetric mode. 
C,,..., C,,together with an additional root at the origin, are the stability roots of 


the aircraft with stick fixed. The appearance of this additional root, and the 
elimination of the root (s— y,), is shown clearly in equations (17) and (18). 


These roots determine modes of motion which are present in all the degrees 
of freedom. The magnitude of the response in each mode, however, depends also 
upon the zeros of the transfer function relating to that degree of freedom. In 
general, any of these roots which is not far removed from an intermediate loop 
pole is not strongly represented in any degree of freedom other than that applying 
to the intermediate loop. For example, the fast subsidence, C,, applies to a motion 
thta is mainly in roll; the oscillatory roots C, and C, apply mainly to yaw, and 
the slow subsidence C, applies to both roll and yaw giving a spiral motion. 


This can be illustrated by considering the response of the aircraft in yaw to an 
applied yawing moment. The closed loop transfer function in this mode is 


[(u, /i.) (s—y)]/Ls (s - (s- « 


ale (s- R,) (s- R,) (s— Rs) (s- Y,) (s- Y,) 


which reduces to 


(42 (s — yv) (s - R,) (s— R,) (s— R,) 
S(s— R,) (s— R,) (s— Rs) (s— Y,) (8 Y.) — (isi) ] x 
x X,) (s- X,) (s— (s- Z,) (s—Z,) 
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The roots of the denominator are C, to C,, together with an additional root at 
the origin, which has come from the transfer function of the forward part of 
the loop. 


The expression can therefore be rewritten 


Hola (s- R,)(s R,) (s — R;) 
igi, — (s C,) (s -C.) (s- C,) (s—C,)’ 


(18) 


in which the term (s—C,) in the denominator has already been cancelled with 
(s—y,) which occurs in the numerator. 


In applying the inverse transformation, the importance of each mode is deter- 
mined by the magnitude of the corresponding coefficient in the partial fraction 
expansion. If any of the roots C, to C, fall close to the intermediate roll poles 
R, to R,, there is approximate cancellation, which results in a near zero coefficient 
for those terms of the partial fraction expansion. This happens with R, and C,, 
so that the rapid subsidence given by C, appears only to a minor extent in the 
yawing motion of the aircraft. 


4.7. AUTOSTABILISATION 


In general, it has been found by experience that the motion which embarrasses 
a human pilot most is the fairly high frequency under-damped mode which is 
primarily dependent on the roots of the simple yaw loop, since the final roots are 
tied to the intermediate roots by the relatively low gain of the outer loop. 


Expression (13) shows that the real part of the simple loop poles is 4 [y, + (n,/i.)], 
and the action of the conventional autostabilising method of artificially increasing 
the value of 7, (i.e. by producing a rudder deflection, and hence a yawing moment 
proportional to rate of yaw) is immediately apparent. 


Reference to Section 4.3, the simple yaw loop (i.e. loop 3 of Fig. 7), indicates 
that this under-damping is due to the smallness of y, and n,/i, and the largeness 
of n,. While it would be undesirable to reduce n, artificially, a beneficial result can 
be obtained by producing a rudder deflection proportional to rate of change of 
sideslip velocity, thus introducing, in the terminology of servo-mechanisms, some 
“lead compensation’’'”. Thus, in Fig. 7 the effect of this autostabiliser is to 
replace n, by (n, + ksn;), where k is the gearing, wherever it occurs. 


The autostabiliser in this case therefore adds a term which increases the 
coefficient of s* in the denominator of the transfer function of the simple yaw loop 
(loop 3), just as with the autostabiliser which operates on yawing velocities. In 
addition, the feed-forward between roll angle and yawing moment is also modified; 
this effect does not seem important, however, and produces a minor reduction in 
the time constant of the spiral divergence. 
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BLOCK DIAGRAM REPRESENTATION 


Ficure 10. The block diagram of a generalised autopilot system. 


It is of interest to note that both these methods of autostabilisation have other 
effects not obvious from this type of analysis. The yawing angular velocity auto- 
stabiliser affects the handling of the aircraft in turning flight, while the rate of 
change of sideslip velocity autostabiliser tends to accentuate the effects of side gusts 
on aircraft heading. 


5. The Aircraft-Autopilot Combination 
5.1. THE GENERALISED AUTOPILOT BLOCK DIAGRAM 


The power of the block diagram approach to aircraft stability problems becomes 
most apparent when it is desired to analyse the behaviour of an aeroplane under 
the control of an autopilot. 


It has been shown how the equations of motion of an aircraft can be represented 
by a block diagram. The additional control loops introduced by a generalised 
autopilot system may readily be added and are shown in Fig. 10, which can be 
simplified to give Fig. 11. Here the autopilot is represented by the blocks A, to A,, 
which represent the shaping networks which modify the aircraft position information 
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Ficure 11. The block diagram of a generalised autopilot system simplified. 


before this is fed to the blocks O, and O,, which in turn convert the autopilot signals 
into aerodynamic moments or forces. Each therefore comprises the transfer 
function of the autopilot amplifier and actuator, and the control surface, or thrust, 
derivative. 


In the lateral case, a conventional aircraft can be controlled only by application 
of rudder and ailerons, which apply moments in yaw and roll respectively. No 
control is present which primarily applies side forces to the aircraft. In the 
longitudinal case control can only be by pitching moment (elevator), and thrust 
(throttles). Thus, in the general lateral case the rudder can be operated in 
dependence upon any function of sideslip velocity, roll angle, and yaw angle, as 
can the ailerons. In the longitudinal case, thrust and elevator angle can be similarly 
controlled by functions of pitch angle, and velocities along, and normal to, the 
aiicraft longitudinal axis. 
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BLOCK DIAGRAM REPRESENTATION 


An aircraft position function measured by any type of instrument can, for 
small changes from steady flight conditions, be expressed in terms of the relevant 
three variables in both longitudinal and lateral cases. Hence, the six blocks A, to 
A,, together with the output blocks O, and O,, when added to the block diagram 
representing the aircraft, describe a perfectly general aircraft-autopilot combination. 
Cross-coupling control derivatives, e.g. i, nz, Y:, Zn, and so on, and the additional 
cross-coupling effects introduced within the automatic pilot, can readily be included 
in the diagram, but have been omitted here for clarity. 


In preliminary studies the diagram can be analysed in stages in a manner 
similar to that already discussed, and the additional terms A, to A, can be selected 
so as to produce the desired modifications to the stability roots of the controlled 
aircraft. In later stages of design, however, the dynamic characteristics of the 
prototype servo-system become available in the form of measured frequency 
responses that cannot be approximated to by simple transfer functions. It is then 
more convenient to turn to the form of analysis based upon the Logarithmic 
Frequency Response, and the Nichols Chart. These techniques are fully explained 
in the standard servo-mechanism literature“. 


Referring again to Fig. 11; the basic control loops | and 3 must first be made 
acceptably stable. Each of these loops comprises the aircraft inertia and damping 
block, and two feedback paths, one aerodynamic and one representing the control 
moments applied by the autopilot in response to changes in aircraft position. To 
reduce a basic loop to a single feedback loop, it is preferable to apply the aero- 
dynamic feedback to the inertia and damping block first, thus producing a stick- 
fixed simple aircraft block (cf. Section 4) upon which the autopilot has to work. 


The stability and performance of the basic autopilot loop can then be investi- 
gated by the Nichols Chart, or equivalent techniques. The autopilot shaping 
networks, blocks A, and A,, and, if necessary, the dynamic characteristics of the 
control mechanisms, blocks O, and O,, may then be adjusted until acceptable 
transfer functions are obtained. 


When satisfactory stability has been obtained in loops 1 and 3, their closed 
loop frequency responses can readily be derived. In addition the frequency 
responses of the feed-forward blocks 2 and 4 may be evaluated by the vector 
addition of the appropriate aircraft and autopilot transfer functions. The product 
of the responses of blocks 1, 2, 3 and 4 is the open loop transfer function of the 
outer loop, the stability of which can be investigated by the conventional methods 
and improved by the adjustment of blocks A, and A,. 


The choice of the remaining autopilot functions A, and A, affects both the 
basic and outer loops. If these terms are required to play an effective part in the 
stabilisation of the basic loops, they have to be accepted in the outer loop. ff 
they are chosen to improve the stability or performance of the outer loop, the 
basic loops have to be recalculated to include their effects. These effects will usually 
not be great. 
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FuiURE 12. Uncompensated response of power control output to a set of low amplitude 
voltage inputs to the autopilot. 


When satisfactory stability has been achieved, a measure of the effectiveness 
of the automatic control may be obtained by plotting the response of the controlled 
aircraft to applied sinusoidal moments or forces, and comparing this with the 
corresponding response for the uncontrolled aircraft. 


5.2. NON-LINEAR EFFECTS 


In practice, the output members of most servo-systems, and in particular the 
hydraulic systems used in powered flying controls, exhibit non-linearities at low 
output amplitudes. The characteristic non-linearity of a hydraulic autopilot 
connected in tandem with a powered flying control system is illustrated in Fig. 12. 
Each full line shows the variation in output phase and amplitude at a fixed frequency 
as the input amplitude is reduced. The set of full lines describe the non-linear 
behaviour over the frequency range. The dotted lines joining points of equal 
input amplitude are the polar gain/phase diagrams at these amplitudes. For a 
linear system the full lines would be radii, since the gain/phase diagrams for each 
input level are identical in form. With the real system, however, the diagram is 
generally distorted so that there is a characteristic “‘swirl’’ near the origin that is a 
measure of the low amplitude phase lag. In the same region a peak is often evident 
on the polar plots, because of lack of damping in the system. 


52 


The Aeronautical Quarterly 


expe 


0: 
LEAD 

toa 

plyi 
thec 

4-0 
| looy 
Ref: 
line 

pos 

At 
of t 
irre 
of 1 

at 
tho 
fre 
inn 
the 

pa 

an 
if t 

6. 
ma 

its 

W 

ch 
the 
dic 

|_| Fe 


erly 


BLOCK DIAGRAM REPRESENTATION 


For consideration of stability, the basic control loops 1 and 3 can be investigated 
experimentally by measuring the frequency responses of the blocks 


to a set of sinusoidal inputs covering the relevant range of amplitudes, and multi- 
plying these by the calculated aircraft response. At each amplitude level linear 
theory can be applied with surprising accuracy to predict limit cycles, stability 
margin at each amplitude, and the closed loop frequency responses of the simple 
loops. Examples of this type of analysis of non-linear elements are to be found in 
Refs. 4, 5, 11 and 12. 


In the outer loop the application of this simple method of dealing with non- 
linearities cannot be applied, since it is not permissible to use the principle of super- 
position in adding into the basic loops the additional signals from the outer loop. 
At frequencies where the loop gain of the simple loop is high, however, the response 
of the aircraft to autopilot input signals will be the inverse of the feedback function, 
irrespective of the non-linearity of the forward path. Thus the transfer function 
of the basic controlled loop will approach 


Ys 
— A3%;) 


Ys 
nd 
at frequencies below their resonance. The linear theory can then be applied to 
those parts of the outer loop due to the autopilot, since in any case the controlled 


frequency of the outer loop cannot be allowed to approach the resonance of the 
inner loop. 


for loop 1, a for loop 3 


The transfer function of the controlled basic loop to aerodynamic moments in 
the non-linear region is more difficult, since the non-linearity is in the feedback 
path in this case. The relative effect of these paths can be calculated for those 
amplitudes at which the control system appears linear, and in the majority of cases, 
if the automatic control is effective, their contribution will be of minor importance. 


6. Conclusions 


An aircraft can be represented by a block diagram, which shows in a simple 
manner the interaction between the forces, moments, and displacements which define 
its flying characteristics. This block diagram can be split into four parts, each of 
which has an obvious physical significance. It has been shown how the dynamic 
characteristics of these elementary parts combine together to give the response of 
the complete aircraft. 


A generalised automatic control system has been superimposed on the aircraft 
diagram, and it has been shown how this can be analysed by conventional servo- 
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mechanism theory and practice. Simple methods of estimating the effects of non- 
linearities in the control system have been proposed, which enable the autopilot 
designer to select the most effective control functions compatible with small 
amplitude stability. 

It is felt that the methods described are of considerable assistance in under- 
standing the full significance of the aircraft equations of motion, and that they should 
be useful both for instructional purposes and in the interpretation of practical results. 
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The Wave Drag of Wing-Quasi-Cylinder 
Combinations at Zero Incidence 


L. E. FRAENKEL 


(Department of Aeronautics, Imperial College) 


SUMMARY: This paper is concerned with the drag, at zero lift, of wing-body 
combinations which are symmetrical about the wing plane. Only two 
assumptions are made: that surface slopes are sufficiently small for the appli- 
cation of linearised theory, and that it is sufficient to satisfy the body boundary 
condition on the surface of a circular cylinder. 


The configuration is represented entirely by singularities along the body 
axis, and three drag formulae are derived. The first involves the Laplace 
transforms of functions representing the strengths of the sources and multi- 
sources; the difficult problem of inverting these transforms is thereby avoided. 
The second involves these strength functions themselves, and is of the form of 
a series of double integrals of the famous von Karman type. The third involves 
functions describing the geometry of a very hypothetical quasi-cylinder which 
has the same drag as the whole configuration in question. The convergence 
of these formulae is established. 


These results are then used to make an order-of-magnitude analysis for 
wing-slender-body combinations: a large number of drag terms are seen to 
be negligible when the ratio of body radius to wing chord is small, and the 
——, terms for smooth bodies are those appearing in Ward’s drag 
ormula“), 


‘ 


Finally, the geometry of the “ equivalent quasi-cylinder ” is investigated. 


Introduction 


A number of methods for calculating the wave drag of supersonic aircraft 
configurations have been proposed recently: Ward’s “area-transfer rule” is 
perhaps the most developed of these. Their great virtue is that the drag is expressed 
as a sum of integrals which involve only the aircraft geometry and certain relatively 
simple “ influence functions.” However, these methods involve two assumptions in 
addition to the fundamental one of small surface slope which permits the use of 
linearised theory. These are that the wing alone or the body alone may be repre- 
sented by a distribution of sources over its surface whose strength is proportional to 
the local streamwise slope; and that the “ interference field,” which must in general 
be introduced to cancel the normal velocity induced by one part of the configuration 
(for example, the wing) at the surface of another part (for example, the body), 
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makes a negligible contribution to the drag. These assumptions can be justified 
in certain rather special cases, but not in general. 


It seems worthwhile therefore to consider a restricted but representative class 
of configurations for which a complete analysis is possible, and to write down the 
“exact ” drag predicted by linearised theory. This is done in the present paper for 
the case of a wing-quasi-cylinder combination at zero incidence: that is, we 
consider a ducted body, not necessarily axially symmetrical, of nearly constant 
radius, and a planar symmetrical wing at zero incidence, the whole configuration 
being symmetrical about the wing plane. 


Three diff‘zent forms of the drag formula are derived: all are algebraically 
(and conceptually) simple but are difficult to evaluate numerically. They are 
intended not to replace the area-transfer rule, but to supplement it. For it seems of 
some interest to display explicitly the terms neglected by this rule: the orders of 
magnitude of these terms, when the body is slender, can then be assessed with some 
confidence. Also, if it should prove possible to evaluate one of these formulae 
numerically, for a few geometrically simple configurations, these results could serve 
as a yardstick with which the accuracy of more approximate methods could be 
measured. Further study of the possibility of numerical calculation is required. 


The mathematical methods of Section 3 of this paper are similar to those used 
by Ward in two unpublished papers which preceded Ref. 1, and which Ward dis- 
carded in favour of that paper. However, the interference field of the present 
analysis was neglected by Ward, and a number of new results are given in the 
present paper. Use is also made of the work of Lance’, Nielsen® and Randall™*. 


NOTATION 
A, (z) multi-source strength function (see equation (15)) 
a mean body radius 
M free stream Mach number 
B=(M? - 1)! 


B,(z) Fourier coefficient of streamwise slope of “equivalent quasi- 
cylinder ” (see equation (16)) 


b wing semi-span, or maximum value of | x | on the wing 
b, (z) Fourier coefficient of body streamwise slope 
D drag force 
G,(z,m) see equation (28) 


*Since the present work was completed, at least two important papers'!!. 12) on the same 
general subject have appeared. However, the present paper still seems sufficiently different 
from Refs. 11 and 12 to warrant publication. 


56 The Aeronautical Quarterly 


« 


le 


ly 


WAVE DRAG OF. WING-BODY COMBINATIONS 


H(z) Heaviside unit function 


r, 


I, modified Bessel function of the first kind“ 

J, Bessel function of the first kind“ 

K, modified Bessel function of the second kind“ 
1 z-value at nose of “ equivalent quasi-cylinder ” 
L_ z-value at base of “ equivalent quasi-cylinder ” 


p_ index of Laplace transform 


S(z) body cross-sectional area 


T ( 


z,x) wing thickness 
z,x) twice the wing streamwise slope, i.e. OT /0z 


U free stream velocity 


V,(z) special function, see equation (22) 


x, 


y,z Cartesian co-ordinates with z-axis parallel to free stream 


Y, Weber’s Bessel function of the second kind“ 


Z(x) z-value on a line of discontinuity of wing slope 


y Euler’s constant 
8, Neumann’s factor: 1 forn=0, 2 forn>1 
¢ measure of surface slope 

n body streamwise slope 

A=F (p) 

w= Bara 

€=BAx 


p, free stream density 


o,(z) see equation (38) 


¢ perturbation velocity potential 


veal and imaginary parts of 


A bar over a symbol denotes a Laplace transform (see equation (2)) 


Suffixes 
b 
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,i,w denote body field, wing field and interference field, respectively 


Nws denotes “ net wing span.” 
57 


8,z cylindrical co-ordinates with z-axis parallel to free stream 
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2. Formulation of the Problem 


Rectangular Cartesian co-ordinates (x, y, z) or cylindrical polar co-ordinates 
(r, #, z) are taken such that the free stream, of velocity U, flows in the direction of 
increasing z, and such that 9=0 on y=0, x >0. The wing has a sharp leading 
edge, is symmetrical about the xz-plane, and has thickness T (z, x) and streamwise 
slope 4t(z, x)=40T/0z. The body has a mean radius a, and its streamwise slope is 


n (z, A= 


bA 


b,, (z) cos né. ; (1) 
0 


The body is taken to be extended from its base to infinity downstream by a cylinder 
parallel to the free stream: sources and multi-sources are required within this 
extension, but their strength decays rapidly as z —> 00. 


The flow field is considered to be the sum of three constituent fields: the wing 
field, denoted by suffix “ w,” which is that of the wing alone; the interference field, 
denoted by suffix “i,” which cancels the normal velocity induced on the body by the 
wing; and the body field, denoted by suffix “b,” which is that of the body alone. 
Each of these fields will be constructed from a distribution of sources and multi- 
sources along the z-axis: the drag of such a line of singularities is then calculated. 
To achieve this, we use the operational calculus in the form given by van der Pol 
and Bremmer?; the definition integral is 


which we also write f (p)=f(z). The modified Bessel functions /, (p) and K, (p) 
play an important part in the operational representation: their inversions are 


1, (-z) H (l-2z*)+A (z- 1), 
0<cos-(-—z) <7, (3) 
(p= sin (n cos~' (- z)) H (1 - z*), 


and K, (p)= sinh (n cosh~'z)H(z—1), cosh-'z>0,_ 


where H is the Heaviside unit function. -(U/2z)f,(p)K,(Bpr)cosn@ is the 
transformed velocity potential of a line of multi-sources of order n, and of strength 
f, (z) per unit length. 


3. The Flow Field 


The wing field is defined by the perturbation velocity potential 


tien t (z,,x,) dz,dx, 


2x) ) /[(z—z,)? - +x,?—2rx, cos 


where the integral is taken over that part of the plan area which lies within the 


(5) 


58 The Aeronautical Quarterly 


fore 


whe 

of . 

apr 

wh 
cire 

anc 

int 
an 

sig 

Fc 

th 

th 

al 

ft 

F 


WAVE DRAG OF WING-BODY COMBINATIONS 


forecone of (r, 9, z), and in operational form this is 
UE. 
Gu 9, P= t (p, x,) K [Bp +x," — 2rx, cos )] dx,, 


where the integral is taken over the span. Let b denote the numerical maximum 
of x, on the wing. The Bessel function K, may be expanded by means of the 
appropriate addition theorem, for r > b, and we obtain 


bw (r, 9, py= - > 5,K,(Bpr) cos n6 t (p, x,) 1,(Bpx,) dx,. (6) 


=0 


where 6, is Neumann’s factor: 5,=1 for n=0, 6,=2 forn2>1. The flow outside a 
circular cylinder containing the wing may thus be represented by a line of sources 
and multi-sources along the axis: the strength of these singularities has a physical 
interpretation as follows. Hayes'’’* has pointed out that, from the viewpoint of 
an observer at (r, 4, z), far from the axis, the sources lying on a “plane of coincident 
signals,” i.e. on a plane 


Bx, cos 6+ By, sin @—z,=constant= - z,, 


may be summed and considered to be concentrated on the z-axis at the point Z,. 
For any given z,, then, these summed source strengths are a function of 6, and 
the strength of the n multi-source in equation (6) is the n™ Fourier coefficient of 
this function. For this multi-source strength is 


| (p, x,) 1,(Bpx,) dx, = | t(p, x,) cos 


o—, 


cos do | t(Z,+ Bx, cos@,x)dx, 


at z=z,. This proves our statement. 


Further, the source strength in equation (6) is Ward’s transferred-area 
function. For, by relations (3) and the product theorem, this source strength is 


Bx |= |2-4| 


which is the function denoted by A’ (z) in Ref. 1. 


*The essence of Hayes’s idea is also reproduced on p. 11 of Ref. 8, which is more readily 
available in the United Kingdom. 
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Consider now the interference field: its potential must satisfy the boundary 
condition 


Ae 2 
6, x,) Kol Bp (a cos . (9) 


where the integral is taken over the wing span. This is the boundary condition on 
the body. The boundary condition on the wing is satisfied automatically; as a 
result of the symmetry imposed, neither the body field nor the interference field 
cause a normal velocity on the wing. So far we have not distinguished between a 
“net ” wing, which is the actual wing outside the body, and a “ gross” wing, which 
is the net wing together with a hypothetical extension of the wing inside the body. 
We now consider a net wing, and denote the “net wing span” by “ Nws”: this is 
the region a= | x | <b for a wing symmetrical port and starboard. In this region 
the Bessel function K, in (9) may be expanded by its addition theorem for | x, | > a. 
and there results 


(a, p)= 2 5, Bp I,/(Bpa) cos né { (en x,)" (p.x,) KABp | x,|)dx, (10) 
NWS 


so that the potential is 


I,’ (Bpa) { 
6 > n 1 
9; (r, 8, p) 5,K,»(Bpr) cos K.’ (Bpa) (sgn x,)"t(p, x,) K,(Bp | x, |) dx 


‘ (11) 
The wing and interference fields therefore combine to the symmetrical form 
U 


1 
by +4 =. K,(Bpr) cos n (Bpa) x 


(sgn x,)"7(p, x,) [ | x, | ) K,(Bp | x, dx, . (12) 


NWS 


The use of a gross wing, of course, does not change this result: extra terms are then 
required in both 4,, and 4;, and they are easily shown to cancel each other. 


The body field must satisfy the boundary condition Ob. =U 5S b,(p) cos nd, 


or n=0 
on r=a, and its potential is 
=U K,(Bpr) cos n Bp K.(Bpa) (13) 
The complete perturbation potential is therefore 
+ G+ > A,(p) K,(Bpr) cos (14) 
n=0 
1 
where A, (p)= (Bpa) 8, (sgn x,)" x,) x 


NWS 


2xb,(p) 
x [ Kx(Bp|x,)) - (15) 
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The functions 22%: Epa) As 


are also of interest. The general behaviour of the functions An (p), Bn (p), An (2), 
B, (z) is not difficult to establish. t¢(z,x) and b,(z) are “cut-off” functions (i.e. 
they vanish outside a finite range of z), and it follows from the definition integral 
(2) that ¢ and b, have no singularities in the finite part of the p-plane; that they are 
O(p) for p—>0; and that for.¥ (p)—> 00 with #(p) bounded, ¢ and 5, are 
O () if t and b, are discontinuous functions of z, O (p~') if t and b, are continuous 
and dt/dz and b,/0z are discontinuous, and so on. It then follows from (15) that 


the singularities of A, are a branch point at the origin, which is also a zero of order 
n+1, and simple poles in # (p)<0. B, has no singularities in the finite part of 
the p-plane; is O (p) for p—> 0; and for./ (p) —> 00 with % (p) bounded, behaves 
like t and b,. 


A, (z) is a multi-source-strength function: it is continuous but vanishes only 
outside a semi-infinite range, decaying rapidly for z —> 00, like z~* for n=0, and 
like z-*"-! for n=>1. (This represents the decay of multi-source strength in the 
cylindrical extension behind the base of the body). B, (z) may be visualised as the 
n® Fourier coefficient of the slope of that quasi-cylinder of radius a which causes 
the same flow as our complete configuration in the region r>b. (The geometry 
of this very hypothetical “equivalent quasi-cylinder” is discussed in Section 6). 
B,, (z) may be discontinuous, but it is a cut-off function. Near a discontinuity of 
B, (z), z=c, say, A, (z) behaves like [z Ba)]}. 


4. The Drag of a Line of Sources and Multi-Sources 


To calculate the drag, we make use of the following form of Parseval’s theorem: 
if f,(z) and f,(z) are two functions with transforms f,(p) and f,(p), whose definition 
integrals exist everywhere on #(p)=0, then 


co ico 

l d 
aD 
— 
provided that both integrals exist. If the integrand on the right has a branch point 
at p=O0, a cut is drawn along the negative real p-axis, and —p is to be taken 
as pe + on the upper and lower halves of the imaginary axis, respectively. 


The drag is calculated from the momentum flux through the curved surface of 
a circular cylinder of finite radius and infinite length, which contains the configura- 
tion. We suppose that all disturbances are confined to z>0. The drag is 
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Deferring questions of convergence for the moment, we substitute 


$= - — A,(p)K,(Bpr)cos nd, 
n=O 
to obtain 


| An(P) Kx (Bpr) Bpr)dp. (19) 
— 

The origin is a branch point and zero of the integrand, and the contributions of the 

upper and lower halves of the path of integration may be combined to give 


D=- 


3 « 1 
poU*Br Ay (pe-*) Eg (Bpr) (Bpre-") + 


n=0 


D= - 


0 


+ K, (Bpre-*) K,’ (Bpr) | dp. 


Now from the relation between K, (p) and K, (pe~), and from the Wronskian of 
I, and K,, it follows that 
K,, (Bpr) K,’ (Bpre-")+K, (Bpre-*) K,’ (Bpr)= - Bor’ 


Also, in view of the definition integral of the Laplace transform, on p=iA 


A, (p) An (pe~**)=| An (id) |. 


2% 
Hence, finally, D= = | A, (iA) | (20) 
n=0 On 
0 
An alternative expression is 


This expression, which is a variant of a formula due to Hayes, can be proved 
equal to (20) by writing the integral as 
co 


0 
iso 


dp 


— iso 


where y is Euler’s constant. The contributions of the upper and lower halves of the 
path of integration may be combined: equation (20) is then recovered. 
A third expression involves the B, (z) and 


K, 
(p) 
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V,,(z) is a function introduced by Ward (Ref. 9, p. 173), and its derivative V,’ (z) 
has been tabulated by Randall. Then 


D=-p,U?xa 


n=0 On 


This expression can be proved equal to (20) by the method already described 
for (21). 


The relative merits of the three drag formulae (20), (21) and (23) depend upon 
the problem at hand. It seems worth emphasising, however, that only (20) avoids 
the vexatious problem of inverting the influence functions which determine the 
A, (z) and B, (z) in terms of the wing and body geometry. 


The convergence of these drag formulae is of some interest, not only from a 
purely mathematical point of view but also from the viewpoint of actual 
computation. This question is discussed in the Appendix. 


5. Wing-Slender-Body Combinations 


For the purposes of this section the wing chord is taken as the fundamental 
unit of length: @ and b are to be taken as dimensionless ratios. It is obvious that 
as the body vanishes (a —> 0), so do the drag contributions of the interference and 
body fields. The rate at which they do so is, however, of some interest, and in this 
section an order-of-magnitude analysis is made for the case of a small. 


We split the functions A, into terms associated with the wing, interference, and 
body fields, and make a change of variable 


2d 
and equation (15) gives, for a half-wing, 
7 
Ais =Od0nl u t ) Jn (€) dé. 
B 
_ (,(%@ & [ _iY la 


6. 
A 
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TABLE I 
ORDERS OF MAGNITUDE OF DRAG TERMS FOR WING-SLENDER-BODY COMBINATIONS 
Dai Div» Dain Dawi 
. Case I—Subsonic wing edges and smooth body 
n>1 ¢,7a* loga € €,,7a7log?a 


Case II—Straight supersonic edges and body not smooth 
n=0 «,? ¢,2a7log?a ¢,¢,a*loga ¢,,2a7log?a 


n>1 ¢,,7a7log?a €,7a? €,,6,@7log a a ¢,,7a7loga 


The order of magnitude of these functions, in terms of » and a, can be 
established in the intervals »<O(a), and the 
process is tedious and will not be reproduced here. The order of the various drag 
terms can then be calculated from (25). 


Table I displays the results of such an analysis. ¢, and ©, denote the order of 
the streamwise slope of the wing and body, respectively, and Dw, for example, 
denotes that part of the drag due to interaction of the body multi-sources of the 
n™ order with the wing multi-sources of the nm order. In Case I the lines of dis- 
continuity of the wing slope t are “ subsonic ” (| dz/dx | > B) and the body slope 
» is continuous; in Case II, ¢ is discontinuous across straight supersonic lines 
(| dz/dx | =contant < B) and 1 is a discontinuous function of z. 


The dominant terms are underlined in Table I: the others are negligible 
compared with these, regardless of the relative magnitudes of ¢y, & and a. If 
‘O (&) ~ O (Ey) ~ O (a), then the dominant terms are O (a*) or O (a* log a). 


The negligible terms in Case I are precisely those neglected by Ward™, who 
also makes the approximation 


A» (2) ~ 27a b, (2) ~ S (2), 
which is justified for smooth slender bodies. 


6. The Equivalent Quasi-Cylinder: the Function Bn (z) 
The function B, (p), which determines the geometry of a quasi-cylinder 
equivalent to our configuration, was defined by (15) and (16) as 


B, (p)= ~ Bp | (sgn x,)"7(p.x,) | (Bp | x, |) Ky’ (Bpa) 


NWS 


(Bpa) Ky (Bp | x, |) 2D 
If G,, (z, m) is defined by 


~ [In (p) Kx’ (mp) - I,’ (mp) K, (p) - () K,, (mp) 1, (mp) K; (p)) 


=G, (p,m), O<m<l, . @ 
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then the product theorem gives 


(sgn x,)" Z-2, 


where f (z,,.x,) vanishes at points off the net wing. Hence, if the interference field 
‘a is to be retained. the usefulness of the drag formula (23) depends upon whether we 


oats are able to invert the influence function G,. We have 

be where & is the Bromwich path c-i°© to c+i0© passing to the right of all 
the singularities of the integrand in the p-plane. 


rag We first deduce the general character of this function. The first term of G, is 
associated with the wing field and the second with the interference field, but the 
+ of properties of G, are simpler than those of either of its constituents. G, has no 
ple, branch point, and for | p | —> © the integrand of (30) is asymptotically 
the 
1 
oa Therefore, by taking large semi-circles in. (p) >0 and & (p) <0, we find that 
ble 
If for z> (l1-m), G,= (m-"+m"), (33) 
the last term coming from the simple pole at p=0. It follows from the definition 
ho integral (2), that the behaviour of f (p) for p—> © indicates the discontinuities of 
f(z); hence we find that 
AG, (- (1—m), m)=AG, ((1 - m), m)= ami’ . (34 


where Af (z)=f(z+)-f(z-). 


The path & may also be replaced by the path ./ , which consists of the 
imaginary p-axis and of an indentation to the right of the origin. The indentation 
ler contributes terms which are half those in (33), and the contributions of the upper 
and lower halves of .V may be combined. Hence 


co 


Guz, m)= sin Az[ J. (md) ¥uQ)] dX (35) 


” This is useful because it gives G, (0, m) immediately, and shows that G, (z, m) - 
G, (0, m) is an odd function of z. 
To invert G, (p,m) completely one may insert integral definitions of the Bessel 
functions into (30) and change the order of integration. The results are seen to be 
8) expressible in elliptic integrals; in particular 
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2 1 1 
G, (z, m)= Jid+m? m (8, k), 


K is the complete elliptic integral of the first kind, and A, is Heuman’s complete 
function’. For n> 0 such results become increasingly complicated. 


The behaviour of the function 


1 Z-Z a 
—— Ga ( 2.) 
| | B\ x, | | x, | 
in the integral of (29) can now be described (see Fig. 1). We suppose z, decreasing 
along a line x,=constant. In region I, z,-—z > B(| x, |—a), the function vanishes. 
In region II, | z,-z|< B(| x, |-a), only the wing term is acting. The function 
rises from 4 (a| x, | )~? on boundary A to 


|2| 
on z,=Z, and rises by an equal amount between z,=z and boundary B. In region 


ill, z,-z<—B(|x,|-a), both wing and interference terms are acting; the 
function is independent of (z, - z) and has twice the value given by (37). 


When z is such that the wing lies entirely in region I or in region III, the wing- 
and-interference part of B, (z) vanishes, (for in the latter case the z,-integration in 
(29) can be made). 


The first term in (37) causes the series 


B,, (z) cos n6é 
n=0 


to diverge, so the equivalent quasi-cylinder, if it exists at all, can certainly not be 


BOUNDARY B BOUNDARY A 


\ 


/ WING 
BODY 


1 / 


FIGURE 1. 
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represented by a Fourier series. But this is not as great a disaster as it might at 
first seem, for by describing the equivalent quasi-cylinder we are merely trying to 


5) attach a physical significance to the function B, (z) appearing in the drag series (23), 
and this series converges. 
te A further insight into the “shape” of the equivalent quasi-cylinder can be 
obtained by examining its geometry at the base. For this purpose we postulate 
that if the wing thickness T (z, x) is not zero at the trailing edge, the wing is to be 
extended to infinity downstream (like the body) by a cylinder parallel to the free 
stream. Let the nose and base of the equivalent quasi-cylinder be at z=/ and z=L, 
respectively, and let 
1g z 
Ta 
5, (z) dz=o, (2), 
(38) 
2 
| By (2) (2) 
yn 
ne Then if Sy is the cross-sectional area of the body of our configuration at the nose, 
Sx+, (z) and Sy +x, (z) are the cross-sectional areas at station z of the body and 
of the equivalent quasi-cylinder respectively. 
“ Now if f (z) is a cut-off function, vanishing outside a< z < b, and if 
@, 
b 
then 7(p)=p F (6)— p* [ BF (b)- { F(2)dz | 
Application of this result to B, (p) gives 
(L)= (sgn x,)" 4 ( | | . G9) 
1 
NWS 


Nw 


Ss 


Equation (39), with n=0, shows that the cross-sectional area of the equivalent 
quasi-cylinder at z=L is that of our complete configuration; equation (40), with 
n=0, shows that its volume is that of our complete configuration. For n> 1, (39) 
describes the distribution of cross-sectional area with respect to @ at z=L. 
Equations (39) and (40) are merely extensions of a result given by Ward” in 
connection with his transferred-area function. 
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7. Discussion 


In this paper an attempt has been made to derive drag formulae which are 
exact within the framework of linearised theory, but which are directly comparable 
with approximate formulae like the area-transfer rule. (In the past, the more 
exact calculations of drag“: '” have largely been based on integration of surface 
pressures, while the approximate ones“':'” have been based on momentum flux 
methods: hence the algebraic formulae have often looked rather different). The 
basic drag formula, equation (20), is a series of positive-definite integrals, each of 
which involves the amplitude of the Fourier transform of a line-source strength or 
line-multi-source strength. This basic formula also appears to be the natural 
mathematical link between drag integrals like (21), which involve multi-source 
strengths and a logarithmic kernel, and integrals like (23), which involve body 
geometry and, in the case of quasi-cylinders, a V,-function kernel. 


It has been confirmed that the difference between the “exact” and approxi- 
mate drag formulae is negligible when the body is smooth and sufficiently slender. 
A largely qualitative description has also been given of the manner in which the 
wing volume would have to be transferred onto the quasi-cylindrical body if such 
an approach were to be used in an exact calculation. 
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Appendix 


CONVERGENCE OF THE DRAG FORMULAE 


General 

If the convergence of the drag formula (20) can be established, the convergence of 
the alternative expressions (21) and (23) is assured, for the individual integrals 
(n=constant) of these expressions certainly converge with respect to z, and are equal to 
the individual integrals of (20). 


We therefore investigate the existence of the double limit 


N->00 n=0 
0 


If this exists the order of taking the single limits may be reversed, because the problem 
is one of absolute convergence. A sufficient condition for the existence of (41) is that 
for A—> oo and/or n—>oo 


For n/A of O (1), that is, in the interior of the first quadrant of the An-plane, A, (iA) is 
therefore required to be O (A~3-«:-®). 


The Wing Field Terms 


The dominant terms in ¢(iA, x) for A—> 00 are due to discontinuities Az; (x) in 
the wing slope across curves z= Z;,(x), say. Such curves of discontinuity will be termed 
“edges,” and we shall consider only the contribution of one typical edge. If ¢ is 
continuous and @t/dz discontinuous, the arguments which follow can be modified by 
replacing Ar by (iA)? Adt/ dz. 


The functions A, (iA) are now split into terms associated with the wing, interference, 
and body fields. Then for a half-wing we have, from (15), 


(in) ~ At (x)exp[-iAZ(x)]J,(BAx)dx. (43) 


n 


2n 


ei"d6 | At (x) exp [- iA{Z (x) - Bx cos 6}] dx (44) 


In the region n>BAb, A, amply satisfies the condition (42); the Bessel function 
J,, is non-oscillatory in this region and 


J,Ay) is O c<cl, for n/y of O(1), r>y, 
ey\" 
and o[n for n>>y. 


In the remainder of the An-plane the behaviour of A,y is closely connected with the 
geometry of the edge z=Z(x). Application of Riemann’s Lemma and of the method 
of stationary phase to the integrals in (43) and (44) yields the following results. 


(i) If the edge is “ subsonic,” | Z’ (x) | > B, the phase of exp[ — iAZ (x)] J, (BAx) 
varies monotonically with x. For A—>00, A,, is O(A~*/?) in general; for 


certain values of n/A, A, is O (A~*/*). 
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(ii) If the edge is “ supersonic” and straight, | Z’ (x) |=k < B. the inner integral 
of (44) is O (A-') for 646, =cos-! (+k/B) and O (1) near 6=8, . (0, identifies 
those planes of coincident signals whose traces on the wing plane are parallel 
to the edge). For A—> 00, Any is O (A~ log A). 


(iii) If the edge is “ sonic * and straight, | Z’(x)|=k=B, then, for A—> 00 


and n bounded, A,,, is O(A~+4); for A—>0o and n/A of O(1), Any is 
O (A? log A). 


(iv) If the edge is “supersonic” and curved, | Z’ (x) | <B, Z” (x)40, then, for k 
A—>oo and n bounded, A,,, is for and n/A of O(1), F 
A, is O(A~') in general. (The inner integral of (44) has a point of stationary 
phase at that x where the trace of the plane of coincident signals associated 
with @ is tangent to the edge). For certain values of m/A, however, 
exp [ — iAZ (x)] and J, (BAx) may have the same frequency as x varies, and 
then A,,, is only O (A-*). 


f~ 


The Interference Field Terms 
In place of (43) we now have 


b 
A nw (id) + An; (iA) ~ 8,,i"-} | At (x) exp [ - iAZ (x)] x 
J, (BAx) Y,’ (BXa) J,’ (BAa) Y, (BAX) 
J,’ (BXa) iY,’ (BXa) 
The Bessel function part of the integrand has the same type of asymptotic behaviour, 


for X and/or n—> 00, as J,(BAx): the preceding analysis therefore covers the new 
terms. 


dx . (45) 


The Body Field Terms 


BX J,/(BXa) iY ,’(BXa) 


If »(z,@) has bounded fluctuation with @, and discontinuities with respect to z, then 
A, is O (n-A-). 

The condition (42) is therefore satisfied in all cases, with the exception of terms 
arising from certain curved supersonic edges at particular values of n/A (see (iv) of this 
Appendix). In these exceptional cases the drag formula can still be shown to converge, 1 
but there are radial lines in the An-plane near which 4A.,,,, is relatively large, and this is ‘ 
undesirable for purposes of computation. 


We now have A,, (iA)= - 4i" (46) 


A Remark on Computation of the Drag Formulae A 
It is clear that in an actual computation the terms associated with the wing field 
would converge much more slowly than the others (see also Section 6). Hence it seems % 
undesirable to use the present drag formulae as they stand: a better approach might wi 
be to calculate, for example, pr 
he 
4x nso 4, | 4. (2)| @)| se 
0 
and, when this reaches a steady value, to add the drag of the wing alone, calculated » 


by some other method. 
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On Damped Free Vibration with Particular 
Reference to Systems having Nearly-Equal 
Natural Frequencies 


R. E. D. BISHOP and D. C. JOHNSON 


(Department of Mechanical Engineering, University College, London; and 


Department of Mechanical Enginering, Leeds University, respectively) 


SumMmMaRY: This paper concerns the free vibration of linear systems which are 
subject to damping. The general theory for a system having n degrees of 
freedom is recalled and illustrated by means of a simple example. Using an 
analytical method, the approximate theory for small damping is derived, 
showing the nature of the entrainment of the principal modes; the results are 
found to confirm those of Southwell (who uses a physical argument) and they 
are illustrated by an example. The approximate theory is shown to need 
modification when there are two nearly-equal natural frequencies and the 
appropriate corrections are found. The analysis is made somewhat com- 
plicated by the existence of two criteria of “ smallness” instead of one; these 
are the difference between the close natural frequencies and the amount of 
the damping. 


1. Theory of Free Oscillations with Friction 
1.1. INTRODUCTION 


This paper is concerned with the effects of viscous damping on free vibration. 
As might be expected, the subject of free vibration with viscous damping has been 
elucidated previously in quite general terms and what is apparently the first complete 
solution of the important problem of small damping has been given by Southwell, 
who uses a physical argument, although an approximate and simple treatment is 
presented by Rayleigh™in his Theory of Sound. Southwell’s results are substantiated 
here by a simple mathematical argument of a quite different sort. A resumé of the 
general theory from which these results are derived occupies the remainder of this 
section. 
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The equations of free motion of an n-degree-of-freedom system which is subject 
to viscous damping are of the type 


Gy Gate. t+ nt + + 
(&=1,2,...H). (1) 


Their solution is given by Lamb’ and it is sufficient merely to summarise the 
results here. 


If solutions qr=V,e™ (3) 
are sought, n algebraic equations are found of the type 


(a, + Dy A V, + + A+ Cop) Vy +. + (AnrA? + Dard (4) 


These relations can be satisfied by non-zero values of the amplitudes WV if 


A,, Aj,...Ajn 
A;, Ax 

A= =0 (5) 


That is to say, A must be a root of an algebraic equation whose order is 2n, and in 
which the coefficients of A are real. 


Corresponding to each root of this equation, there is a definite ratio between 
the W’s. This ratio is conveniently expressed in the form 


WV, 
(7) 
a, An 

where 2,, 2... . 2, are the cofactors of any one row in A when A takes the value 


of the root concerned. 


In examining this result, attention must be paid to the form of the roots. They 
may occur in pairs (A,, A’,), (A,, A’)... - (A,,, A’,), the members of each pair being 
complex conjugates, or they are real and different*. 


tSee also Duncan), 


*The question of repeated roots is of mathematical interest only and will not be raised here. 
The more important matter of close roots will be mentioned later. 
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ot Thus there may be an s“ pair of roots which can be written as 

A,=pstio, 
l) @) 
») Alternatively the roots may be real numbers, the s pair of roots being 
) The identity of these real roots as a pair can be maintained, since a pair of the 


form (8) disappears when two real roots appear, as the damping is increased 

proportionally throughout the system. Thus, of the 2” roots of equation (5), a 

number may exist of the type (8) (in conjugate pairs) and the remaining pairs of 
) roots will be real of the form (9). 


These two types, which correspond to the conditions of “‘light’” and “‘heavy” 
(or “over critical”) damping respectively, will be discussed separately. 


NOTATION 
A,, see equation (6)* 
d,s inertia coefficient for generalised co-ordinates q* 
a, inertia coefficient for principal co-ordinates p* 
) b damping coefficient (see Fig. 1 on p. 78) 


b,, damping coefficient for co-ordinates p and gq; the same symbol is 
used for both sets of co-ordinates, although the actual values (for a 
given system) will be different* 


’ :C, modulus of cofactor defined in equation (10)* 
C,; Stability coefficient for generalised co-ordinates g* 
c, Stability coefficient for principal co-ordinates p* 
G,,H, constants of integration for st mode* 
h_ see equation (70) 
h,,h, roots of equation (72) 
I polar moment of inertia of disc (see Fig. 1 on p. 78) 
K,k torsional stiffnesses of light shafts (see Fig. 1 on p. 78) 


n number of degrees of freedom 
Pp, principal co-ordinate* 


February 1958 


R. E. D. BISHOP AND D. C. JOHNSON 


qr generalised co-ordinate* 
R, R’ see equations (34) and (36) 


a2 dimensionless measure of damping defined in equation (28) for the 
system of Fig. 1 (p. 78) 


2, cofactor of r element of any row of A withA=A, 
B,, Ys s™ pair of real roots of equation (5) 
IT, I” see equations (34) and (36) 
A see equations (5) and (17) 
€=C,-—C, for a system in which a,=a,; see equation (68) 
6, argument of cofactor defined in equation (10)* 
A see equation (30) 
A see equations (3) and (15) 
A,, A’, complex conjugate roots of equation A =0; see equation (8)* 
II, amplitude of p,* 
Ps,o, see equation (8)* 
Ww, amplitude of g,* 


©, natural frequency* 


1.2. LicHt DAMPING 


The pair of roots (8) give values of the cofactors a,, 2, ... %, of a particular 
row of A which are themselves complex conjugates. The notation will thus be 
adopted: — 


apm fC, (for A ==A,) (10) 
(for J 

If now the pairs of results (8) and (10) are taken back to the trial solutions (3), the 
latter are found (after rearrangement) to yield an s‘* mode of oscillation of the type 


cos (c,t + +H, sin + 9,)] 


La” [G, cos (ot + +H, sin + .9,)] 
(11) 


This type of motion will be referred to hereafter as that of a “damped mode.” 
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DAMPED FREE VIBRATION 


A question of nomenclature arises here. Undamped systems are frequently 
said to possess “normal modes” of vibration. These modes are defined for a 
system, then, in the absence of its damping. The authors prefer to use the adjective 
“principal’’ which is also commonly used to describe these modes of conservative 
systems. The word “normal” can then be reserved to indicate that a co-ordinate 
has been selected in a particular way to describe the intensity of distortion in one 
of these modes. 


Now a damped mode is unlike a principal mode of oscillation in that there 
are phase differences within it, as shown by the constants .6,, ,0,, and so on. 
Moreover, the amplitude diminishes exponentially owing to the presence of the 
factor e”*’ in each expression; a proof that p, must be negative for a stable system 
is given by Lamb’. This mode is evidently independent of the remaining damped 
modes and motion in it decays independently of motion in the others. 


The distinction between a damped mode and a principal mode is a fundamental 
one and it is this which will be examined in subsequent sections. There is one 
difference, however, which may be noted immediately. Owing to the phase 
differences within any one damped mode, no fixed distortion can be regarded as 
being made up of one damped mode only. It is therefore not possible, in general, 
to start a vibration in one damped mode only by releasing the system from some 
particular static distortion. 


It should be noted that, strictly speaking, the words “‘frequency”’ and “‘phase” 
should not be applied to the quantities « and 4, since the motion is no longer 
harmonic. However, these terms will be retained, as confusion is unlikely to arise 
and suitable substitutes are not easy to find. 


If all the damping coefficients b,, are reduced to zero, then it may easily be 
shown that, in equations (11), 


and that the constants .C, tend to values which determine the s‘" principal mode. 


1.3. DAMPING 


If the damping forces are sufficiently great within a system, they may mask 
those of inertia and elasticity and the motion may become aperiodic. There is then 
an even number of roots of the type (9) and each yields real cofactors z, which may 
be written 


(for A=A,) 


(for A= r’,) 


February 1958 75 


R. E D. BISHOP AND D. C. JOHNSON 


The results (9) and (13), when taken back to the original trial solutions (3), yield 
an s mode of the type 


q,=G, +H, .E,e™ 
q.=G, +H, .E.e™ 


B 
qn=G,; sD,e + H, 


where the constants G, and H, may take any real value. As might be expected, 
8, and y, are negative for a stable system. 


It will be seen that these last two components could be treated as separate 
“‘modes’’; but this would involve a change in the total number of modes as the 
damping is increased and hence some of the simplicity of the treatment would be lost. 


Generally speaking, this type of motion is not of great interest to the vibration 
analyst as it is observed only in systems with sufficient damping to make violent 
oscillation unlikely. 


1.4. SUMMARY OF THE FOREGOING RESULTS 


To sum up then, damped free vibration occurs in modes which may either be 
of the oscillatory type (11) or of the aperiodic type (14). If all the damping forces 
are light, then all the modes are oscillatory; but if the damping is (unusually) great 
some, at least, of the modes will be aperiodic. When the motions in these modes 
are added together, there are 2m constants (which have been denoted by G, and H,) 
and these are sufficient in number to allow initial conditions on the qg’s and q’s 
to be fitted. 


If the system oscillates freely in the mode (11) and the damping is fairly small, 
then although motion clearly does not take place in any single principal mode, it 
is to be expected on physical grounds that the s“ will predominate. Coupling—or 
“‘entrainment”—occurs between the principal modes. The nature of this coupling 
appears to have been discussed first by Rayleigh’), whose approximate treatment 
forms the basis of present engineering practice. Using a simple argument, Rayleigh 
made certain observations whose correctness has been confirmed by Southwell”’, 
who gives a more complete treatment which is based (like Rayleigh’s) upon a physical 
argument. The alternative treatment of Section 4 arrives at Southwell’s results by 
simple mathematical considerations. 


2. The Use of Principal Co-ordinates 


Since the principal co-ordinates of a system are merely a special set of 
generalised co-ordinates, the foregoing theory of Section 1 still holds good when 
they are used. Two features of the principal co-ordinates make it worth while to 
adopt them before proceeding. In the first place, the coefficients a,, and c,, (which 
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correspond to them) all vanish (r ~ s) although the coefficients b,, remain finite, in 
general (but with different values from those of any other set); thus principal 
co-ordinates p,, Pp. . ... Pn reduce the amount of algebraic manipulation which 
must be performed. Secondly, the value of each principal co-ordinate is a direct 
measure of the intensity with which the appropriate principal mode is present in a 
given distortion*. This paper is thus concerned with the relative magnitudes of the 
co-ordinates p during motion in a damped mode. The same analysis as before is 
used at first, but with the co-ordinates p. 


If trial solutions . ‘ (15) 


are introduced into the equations for the principal co-ordinates, n algebraic relations 
are found of the type 


b, AIL, +...+(a,A? +b, (16) 
That is to say, it is now necessary that 


A, 


where A,=aA? + b,.A+C,. ‘ ; (18) 


There are 2n roots of this equation which are the same as those found previously 
and for light damping these occur in conjugate pairs. Corresponding to any one 
root, there exists a set of ratios between the amplitudes II, which are given by 


(19) 


where z,, 2... 2, are the cofactors of any one row in A when A takes the value 
of the root. 


3. Example 


The results of Section 2 may be illustrated by reference to the simple torsional 
system shown in Fig. 1, in which rigid flywheels perform torsional oscillations with 
torsion springs of stiffnesses k and K providing restoring couples. Let the displace- 
ments of the discs from the equilibrium condition be qg, and g, as shown. The 
system is damped at g, by a viscous damper whose damping coefficient is b. 


*The scale of measurement which is used is a matter of choice, so that the intensity of distortion 
in the r* principal mode, which is meant by the statement that p.=1, is left open. Normalised 
principal co-ordinates are a set of co-ordinates for which this choice is made in a particular 
way. 
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FiGure 1. 


If there were no damping present, the system would be capable of oscillation 


in the two principal modes: — 


Thus principal co-ordinates p, and p, may be taken such that 


Pe 
With this set, we have 
+ 


In terms of the co-ordinates p, the equations of free damped motion are 


21, +bp,+2kp,—bp.=0 
—bp, +2 


The trial solutions 


Pi 
p,=IT,e* 


now give a pair of homogeneous algebraic equations, namely 


E + (2) II, - (5) AIT, =0 


-(%)an, + (F)a+o, | II,=0 
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(20) 


(21) 


(22) 


(23) 


(24) 


(25) 
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DAMPED FREB VIBRATION 
where the quantities , and , are the natural frequencies of the system in the 
absence of its damper. That is 


k+2K 


(26) 


To illustrate this theory, we now trace the effects on the solutions (24) of 
increasing the damping and, for definiteness, we shall take 


K=2k, so that o,?=Sw,?.. ‘ (27) 
n 
The computations are greatly simplified by introducing a dimensionless parameter 
2 as a measure of the damping such that 
)) b 
These values may be substituted into equations (25) and the determinantal relation 
corresponding to equation (17) written down; this gives a quartic in A,. namely 
42a0,A? + 60,202 . . . (29) 
This equation can now be made non-dimensional by the substitution 
A 
A= ; ‘ (30) 
) 
which yields . (31) 
The roots of this equation can be found for any value of z. By substituting roots 
back into either of equations (25), the corresponding ratio of the amplitudes II, 
and II, may be found; the variation of the roots of the equation has been plotted 
) in Fig. 2. 
Consider first the roots of equation (31). When 2=0, they take the values 
A=i, -i, -i/5, . ‘ ‘ (32) 
) in conformity with the simple theory for the conservative system; these pairs of 


roots are represented by points A and G in Fig. 2. 


If « is increased, the roots form two pairs of the type (8). Their imaginary 
parts become smaller than the values (32) and they are represented by the drooping 
curves AB and GH, which refer to the dimensionless quantities o,/, and o,/,. 
As « is increased from zero, the roots A acquire real parts p,/, and p,/, which are 
negative; these quantities are represented (in the positive sense for convenience) by 
the curves OC and OF respectively. 
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FIGURE 2. 


When a= 1-474, the pair of complex conjugate roots -\,, \’, (which corresponds 
to the first damped mode) changes over to form two real roots of the type (9). 
These are negative and are given by the curves CD and CE, which represent the 
dimensionless quantities - 8,/, and —y,/,. This system of roots—that is, two 
complex and two real—persists for all values of z greater than 1-474. 


If both flywheels had been provided with dampers and the two damping 
coefficients were both increased, all motion would become aperiodic so that 
ultimately all four roots of the appropriate equation in would become real. This 
is not so here, because an unlimited increase of 2 eventually fixes the right hand 
disc, leaving a conservative system with only the freedom g,. When z is large the 
second and fourth terms in equation (31) predominate; and, if the others are 
discarded, the root +i/3 is found immediately. That is to say, when z—> 00, 
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iV 
(33) 


This is clearly correct, as the left hand disc / of Fig. 1 is left to oscillate freely, the 
total stiffness of the controlling springs being 3k. Moreover, as 2—>0o, both of 
equations (25) show that II, and II, approach equality so that, by equation (21), 
q, diminishes indefinitely whereas q, does not. 


The relative amplitudes and the phase difference between the distortions in the 
principal modes during vibration in a single damped mode may be found from 
equation (25). They are R and I’, respectively, where 
Il, A?+eA+1 
= Re’, ‘ (34) 
A taking the value of the appropriate root of equation (31)*. R and I have been 
evaluated for different values of 2 using the root 


(given by the curves OC and AB of Fig. 2). The variation of R with 2 is shown 
in Fig. 3 and it will be seen that, as might be expected, R=0 when 2=0; this 
curve is terminated at «= 1-474, since the first damped mode ceases to be periodic 
at that intensity of damping. The value of the phase angle I is also plotted in 
Fig. 3 and is seen to vary between 7/2 and =. 


It can easily be shown, in fact, that when damping is increased within any 
system to the point at which a damped mode becomes aperiodic, then the motions 
in the oscillatory principal modes just before they vanish (which together constitute 
the motion in the damped mode) are either in phase or in antiphase with each other. 


The ratio of the complex amplitudes II for the second damped mode is more 
conveniently found in the form 


Il, _ A?+aA+5 


— Rail’ 
il, = R’e™ (say) (36) 


where A now takes the value 


*The two members of a pair of complex conjugate roots produce complex conjugate values for 
this ratio in the form Re*'. When the result is used in equations (24) and the changeover 
is made to trigonometric functions in the form (11), I is found as a phase difference between 
p, and 
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FIGURE 3. 


Using the positive sign in this expression, the variations of R’ and IY with z have 
also been computed. They are also represented by curves in Fig. 3. Once more 
the curve of the modulus of the ratio (i.e. of R’) passes through the origin, as is 
demanded by the simple theory of conservative systems. The angle IY varies from 
— 7/2 to zero and it is seen that the complex ratio II,/II, tends to the value unity 
as 2 tends to infinity, in accordance with the previous prediction. 


Before leaving this example, it is worth pointing out a feature which is of 
practical significance. It is that, when the damping is not small, its effects upon 
the motion are both considerable and complicated. It would generally be difficult 
to deduce the principal modes of a fairly heavily damped system from observation 
of a free vibration, even if the latter occurred in a single damped mode. It would 
also be difficult, in general, to predict the variation of the entrainment of the principal 
modes within a damped mode as the intensity of one or more damping forces is 
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varied. As an example, it might be mentioned that these considerations should be 
borne in mind when free vibrations are observed of airframes which are supported 
on partially deflated tyres or on air-bags on the ground. 

4. Small Damping 


It will be found that the determinant A of equation (17) is of the form 
M=A,A,...An+O[b?, A], (38) 


where O [b?, A] signifies terms which contain the cross-coefficients b,,.(uv) to the 
second power at least and are also functions of A. 


The roots of equation (17) are evidently dependent on, and therefore functions 
of, the coefficients b. Equation (38) shows that if only the coefficients 5,,, 5... .. Dan 
were present, the coefficients b,,, being zero, then the roots would be 


If the direct coefficients b are small, so that their second powers may be neglected, 
these roots may be written approximately, as usual for light damping, 


(39) 


b b 

A,=-x' A’,=-  -io 

1 2a, 1 1 2a, 1 

b 

Ban Dia 

2a), n 

where w,, , . . . , are the natural frequencies of the undamped system. 


It is now necessary to examine how these roots will be modified by the presence 
of the cross-coefficients b,,. To do this, consider a particular root A, and let it 
be expressed in ascending powers of the coefficients b,,, using Taylor’s theorem. 


Evidently A,= - 3 [ Bue | +... (higher degree terms), (41) 
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and, when this expression is evaluated, the partiai derivatives must all be related 
to the condition where 


bw=0 (ux) 
(42) 


8 


+i, 


Now A, is not given as an explicit function of the coefficients b,, and the 
differentiation cannot, therefore, be performed directly. Instead, the relation 


(uAv). . . (43) 


must be used, because the function A relates the b’s and the root A,. But it has 
already been pointed out that in equation (38) the b,, terms in the expansion of A 
are all of the second degree or higher. When A is differentiated with respect to the 
coefficients b,,, therefore, these terms all vanish if b,, is put equal to zero. The 
numerator of equation (43) is therefore zero. The addition of cross-damping terms 
to the equations thus produces second-order changes only in the values of the 
roots of equation (17) if the coefficients b,,, are of the first order of smallness. 


Any pair of roots A,, A’,, as given by equations (40), may now be substituted 
into the equations (16) in order to find the time-variation of the damped motion. 
During this motion the ratios of the displacements at the different co-ordinates is 
given by the co-factors of A, as in equations (19); and, in finding these ratios for 
the root A, it is convenient to use the co-factors of the s‘" row. The co-factor of the 
‘*“ term of this row will now be found; this is the quantity z, of equations (19). 


When finding the co-factors of A it is necessary to write A, for A. With this 
substitution it can be seen, from an examination of A, that 2, can be written in 
the form 


a, =(@,A,7 + B, +C,) + + Ce). (AnAg? + DanAs + Cn) ox +O [b?, Ag] (44) 


where the suffix “‘Os”’ after the product indicates that the term (a,A,* + b,.A,+C,) is 
missing, and where O[b?,A,] denotes terms containing A, which are of the second 
degree or higher in the coefficients 5... If the b,, terms are sufficiently small, this 
last expression O [b, A,] will be negligible, so that only the product is left. This, 
then, with A, made equal to 

b 


a, +i0,, . : ; (45) 


gives the required cofactor 2, for use in equations (19). Similar treatment holds 
good for the conjugate root \’, and it will be seen to give the conjugate result for 2.. 
84 
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The other co-factors of the s row of A are now required. Consider the r® 
of these, which has been denoted by z,. This determinant can be expressed in the 
form 


where 

eee bends | 
bo A; ba As (a,A,? + + C2) 
| 

(47) 


This arrangement of the determinant is obtained from A by moving the s column 
to the first column and the r” row to the first row. This gives b,,A, as the first 
term of the first row, and there is no other 5,,A, term, because it is removed in 
forming the co-factor from 4. It can now be seen that 2, can be put in the form 


tp = By As [(Q,Ag? + By +C,) (AAs? + + DanAs + Cn) Joros + OLB*, Ac], 
(48) 


where the subscript “‘OrOs” indicates that the factors (a,A,?+5,A,+c,) and 
(a,A,? + b,,A, +¢,) are missing from the product and where O[b*, A,] denotes a function 
of A, whose terms are of the third degree or higher in the coefficients b,,. If the 5,. 
are sufficiently small, this function can be neglected in comparison with the product 
function, although it should be noted that this is itself a quantity of the first order 
of smallness, since it contains the factor 5,,A,. 


It has thus been found that 


@,=(a,A,? +B, +) (GAs? + t+ Co). (AnAs? + DanArs + Cndos (49) 
A, [(a,A,’ b, iAs +¢,) (a,A," C2) ees (a,A,” + Bands + Ca)Joros 


These results correspond to the root 


of equation (17) and are only appropriate when the coefficients 5 are sufficiently 
small. If the conjugate root 
_ Oss 

2a, 


-iw, . ‘ ‘ ‘ (51) 


is used, then the quantity A, in equations (49) must be replaced by 2’.. 
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The variations of the principal co-ordinates p, and p, can now be written down 
by the use of equations (19). If the ratio II,/II, is formed, it is found that 


a, b, A, 


When the value of A, is inserted from equation (50), and terms of higher degree 
than the first in b are neglected, this becomes 


Tl, 


The ratio is thus pure imaginary, and if \’, is used instead of A, the conjugate 
result is arrived at. And when the results corresponding to A, and A’, are taken 
together, in the manner of equations (11), it is found that, if 


p.=e 1G, cos w,t +H, sin o,f], | (54) 
: 
then [G, cos (w,f +42) +H,sin(w,t+4)] . (55) 


The factor b,, in the numerator of equation (53) shows that the fraction is 
small and therefore the amplitude of the r“ principal co-ordinate is small compared 
with that of the s‘ principal co-ordinate when the system is vibrating in the s" 
damped mode—that is when the motion corresponds to the s“ pair of roots, A, 
and ’,, of equation (17). In the analysis p, is any principal co-ordinate other than 
the s‘, and the result shows therefore that all the principal co-ordinates are small 
compared with the s‘“". And, in general, this effect becomes more and more marked 
as the remoteness of », from w, increases, owing to the appearance of the factor 
c,-a,,’; for c,/a,=,*. Further, all these co-ordinates p, are in quadrature 
with the s‘, since the ratio II,/II, is imaginary. 


These results may be illustrated by reference to the system of Fig. 1. The 
curves of o,/, and +,/, shown in Fig. 2 have zero slope when z is zero, so that. 
for small damping, the frequencies +, and o, are very nearly equal to ihe natural 
frequencies ©, and »,. The curves OC and OF of Fig. 2, on the other hand, have 
a negative finite slope at the origin, so that, to the first order, a decaying motion 
is predicted; these two curves happen to have the same slope at 2=0 because 
b,,=,,. Figure 3 shows that the entrainment of the principal modes is small 
when the damping is small and that the entrained principal mode is nearly in 
quadrature with the predominating principal mode in either damped mode. 


It is now necessary to return to the question of how large the coefficients 5 
may be, in comparison with the other energy coefficients, before these approxi- 
mations cease to be of value. Firstly, it is necessary for the s‘ mode to be lightly 
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damped; otherwise the approximation for A, in equation (40) breaks down. The 
approximation of that equation requires that 


Secondly, the cross-coefficients b,, must be sufficiently small for the second-degree 
terms in equation (41) to be negligible in comparison with b,,/2a,. In discussing 
this condition it is convenient to assume that the co-ordinates are defined in such a 
way that all the inertia coefficients a are equal. The condition will then hold 
provided that any b,,.A, term in equations (16) is small compared with a@d,?._ That 
is to say, if 


5 


for all the cross-coefficients b,,. 


Other conditions must be fulfilled if the relation (53) is to be valid. So far no 
restriction has been placed on the magnitude of the coefficients b,,; if one of these 
is large, however, so that b,, is not much smaller than a,o,, then the factor ,, 
which contains it, will contain an appreciable imaginary component and so p, and 
p, will not be in quadrature. 


These conclusions can be expressed in another way. If all the inertia coefficients 
are equal to a, then the damping may certainly be treated as “‘small’’ if 


b | b 


where b represents any damping coefficient, either direct or cross. 


Finally, a new condition is introduced by the denominator of equation (53). 
If , is very close to ,, then this denominator will itself be small and the amplitude 
of the r co-ordinate p, may then be comparable with that of the s‘" for motion in 
the s‘ damped mode. This condition may also interfere with the phase relations 
and it is discussed in the next section. 


This deduction of the form of the free motion in a lightly damped system may 
also be made by a simple algebraic argument as follows, starting from equation (16). 
The damping is supposed to be light, so that the coefficients b are small. Now 
consider solutions which only differ from those which are valid when the b’s are 
all zero by having in any mode, say the s, small components of .. . 
in addition to a large II, component. These other small components can be 
neglected, in comparison with II,, in the s‘® equation; this then becomes 


(a,A? + b,A ,=0 . (59) 


and determines the two roots of A, as previously discussed. If one of these roots is 
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now substituted into one of the other equations, say the r, and if in this equation 
all the b terms are neglected except that containing IT,, then it is found that 


. . . (60) 


° II, br Xs 


This equation corresponds to the relation (53) and illustrates the quadrature 
phase difference between the II, and Il, components; for the numerator b,,, has a 
large imaginary part compared with its real part, whereas the denominator is 
mainly real. This general argument was given by Rayleigh”. It reveals very 
simply how the relationships between the II’s come about, but it is not a guide 
for the limits of application of the approximate theory. The foregoing theory is 
more complete and provides this; and, as shown in the next section, the present 
theory can be extended to the case of close natural frequencies. 


It is also possible to explain the behaviour of the system by direct argument as 
follows. If the coefficients b,,(r 4s) were absent, then the system could vibrate 
freely in, say, the s‘* mode with the appropriate frequency and damping factor. 
The addition of a single damping element, which would introduce a small b,, term, 
would cause a force to be transmitted to the r’ mode because of the motion in the 
s". This force would be in phase with the rate of change of the displacement causing 
it, and the displacement which it would itself cause would be in phase—or anti-phase 
—with the force, provided that the r frequency were not very close to the s and 
the r‘" mode were not very heavily damped. The quadrature r* component thus 
introduced would not appreciably affect the s“* frequency or damping, provided 
that 5,, is sufficiently small. 


5. Nearly-Equal Natural Frequencies 


The conclusions of the preceding section are invalid if two natural frequencies 
of the undamped system are very close. This is because the denominator of 
equation (53) may then become so small that it is comparable in value with the 
numerator. Thus if , is very close to w,, the amplitude of the r principal 
co-ordinate may be comparable with that of the s co-ordinate, when the system 
is vibrating in the s‘* damped mode. This result is not surprising, since motion at 
the r co-ordinate p, is excited, through the damping forces, at a frequency which 
is close to resonance, so that the amplitude may be large. The proximity of the 
r™ and s‘" resonance peaks also suggests that the phase relations between the r" 
and s* modes will not be the same as they would be if , and , were not close to 
each other. It will be shown that this may be so and the amplitude relations will be 
examined more fully. Finally, it will be demonstrated that small changes in the 
frequency of free vibration are produced by the damping when two of the undamped 
frequencies are very close, the changes being such that the frequencies are drawn 
closer together; this will be done by treating a simplified case. 
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Consider the motion in the s“ mode, when the frequency w, is very close to «,, 
and let the damping of the r' and s‘" modes be small, so that equations (40) are 
still valid when the cross-coefficients 5,,, are zero. For the present it is convenient 
to continue to assume that, even when the coefficients 5,, are not zero, all the 
second-degree terms in the b’s, in the expansion (41) of A,, can be neglected; but 
it will be necessary to return to this matter later. 


The expressions for 2, and a, determine the relative phase and amplitude ratio 
of the r and s‘ principal co-ordinates during motion in the s‘* damped mode. 
The expressions are those of equations (49) when 


8s 


+10, 


is substituted for A,. The propinquity of w, and , affects the value of the factor 


2 


Dz, ° bes at ) 
a, ( 2a, + b,, 2a, J (62) 


which occurs in 2, but not in z,. This factor also occurs in the expressions for 
G1, %... Sy, Z4,, and so on. If second-degree terms in the b’s are neglected, 
the factor becomes 


a, 


(C, — +1 bss Js, (63) 


from which it can be seen that the real part may not be very much greater than the 
imaginary part if c,/a,, that is to say, ,” is very close to w,’. 


Consider the extreme case when 


(c,-—a,02) << - bu.) Os, (64) 


so that the expression (63) is almost a pure imaginary. In this event the result (53) 
may be replaced by 


i, bys 
& 
b,,- = bes 


where p, is any principal co-ordinate other than the r” and s‘*. Thus the amplitudes 
of the r and s“ co-ordinates are comparable and the corresponding distortions 
are in phase. The amplitude of any other co-ordinate is small compared with them 
and the corresponding distortions are in quadrature. 
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The quantity 


which occurs in the denominator of the expression for II,/IT,, will be zero if 


This relation is likely to hold good if the two modes whose frequencies are close 
are similar in shape, and this is usually the case in real systems where close 
frequencies are met with. It applies, for instance, to straight circular bars which 
have flexural modes in two planes at right angles. For this type of system, the 
assumed inequality (64) is not likely to hold and the nature of the motion must be 
examined further. 


A full investigation of the free motions of a system with close natural frequencies 
cannot be made without allowing for the effects of the second-order terms which 
were previously neglected in the Taylor expansions. To attempt to introduce these 
terms into the general theory would involve excessively complicated analysis. This 
is because there are now two criteria of “‘smallness” to consider instead of one, 
namely the difference between , and w, and also the amount of the damping; 
further, the introduction of specific second-order terms into the Taylor expansion 
leads to excessive complication. Instead of doing this, therefore, it is profitable 
to examine the theory for a lightly damped system with two degrees of freedom 
only; this will illustrate the main features of the motion of more complex systems. 


6. Example 


Let the principal co-ordinates of a two-degree-of-freedom system be p, and p., 
so that the equations of motion are 
py +Cip, +b,.p,=0 
. (67) 
Py +2 P2t+ bog 


The co-ordinates may be chosen in such a way that a,=a, (as with the system of 
Fig. 1). Let the stability coefficient c, be equal to c,+¢ where ¢/c, is small 
compared with unity, so that », is nearly equal to w,. Suppose further that the 
damping is small; that is to say, 


where is the absolute value of b,,, 6.. or 
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The determinant A for this system will have two conjugate pairs of roots and 
both pairs will be approximately equal to 


b 
+10,, 


which represents the approximate pair of roots for the first damped mode when 
@, is not nearly equal to w,. Taking the positive sign, a solution may be sought 
of the determinantal equation 


lar? +b,,A+c, by | 
| =0, (69) 
b,,A + | 
which is of the form h=A,;+h (say) ‘ 4 (70) 
1 


where h is real or complex, but has components both of which are small compared 
with »,. If this trial solution is substituted into equation (69), and all terms of 
higher degree than the second in 


a,®, 


are neglected, then it is found that 


This equation is a quadratic in h/w, and it has the roots 


The roots, h, and h, say, are complex in general: The roots for the free motion 
of the complete system therefore differ in both frequency and damping factor from 


_ those which would hold for the first principal mode if it were entirely independent of 


the second. For both the real and imaginary parts of equation (70) contain 
contributions from h. 


The modal shape for the frequency A, +4, can be obtained from the relation 


— A, +h,) 
2 4,(A,+h,)?+5,, A, +h) +e, 


+h,)? + (A, +0¢,+¢ 
—b,, (A, +h,) 
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and similarly for A,+h,. And, as both the real and imaginary parts of A are 
small compared with the imaginary part of A, this relation can be written 
approximately as 


II,” 2h, /0, b,./(a, y 


An idea of the form of the motion in particular cases can be obtained from this 
expansion. 


Before examining these results, however, it is worth going back to equation (69) 
to seek a solution of the remaining form 


The small quantity h’ may now be evaluated, just as before, and, in fact, is governed 
by equation (71) with i replaced by —i. The roots h’, and h’, are thus the complex 
conjugates of h, and h,. If A, +h, is replaced by A’, +h’, in equation (73), then 
the conjugates of the previous expressions for II, /II, are arrived at. The frequencies 
and amplitude ratios are therefore unchanged and only the phase relations are 
altered. 


Consider first the case when b,,=b,, which, as mentioned earlier, is common 
in practice. Equation (72) now reduces to 


4b,. 


4 a,?0,? 
bi, 
In this event, equation (76) gives 
Ah ie 5b,’ ©, 
Substitution of the first of these values into equations (74) shows that 
 _. 
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where R >> 1, while the second gives 


II, b,2/(a,,) = i 
I, e/c, R° 


(80) 


Now the meaning of these results can best be seen when it is remembered that, if 
b,. is zero, then the roots of equation (69) are 


jo, (1+ = =): 


The two imaginary values of A which have just been found for small 5,, show that 
the natural frequencies differ from these. Thus the lower frequency is raised by 


_ 


and the higher frequency is lowered by the same amount; the 5,, factor has the 
effect, therefore, of drawing the two frequencies very slightly closer together. The 
expressions (79) and (80) show that the two principal modes are almost in quad- 
rature during motion in either damped mode; at the lower frequency, the amplitude 
of the first principal mode is much greater than that of the second, whereas the 
opposite holds at the higher frequency. The motion in the second principal mode 
lags behind that in the first. 


If the calculations are repeated, using the conjugate root A’, in place of A,, 
then the frequency and amplitude relations are unchanged, but the phase relation 
alters from lagging to leading. 


The motion of the system can be depicted as that of a particle moving in a 
plane. The two free motions of the particle, when b,, is zero, correspond to two 
motions along perpendicular directions. When 6,, is introduced, the motions 
become ellipses which diminish in size as the motion is damped out. The ratio of 
the axes of the ellipses is given by 


e/c, 


b,,/(a,,) 


If ¢ is very large compared to (w,b,,), then the ellipses are long and thin; they 
become progressively rounder as < is reduced. If < is made too small, then the 
approximations in the foregoing treatment become invalid and the results do not 


apply. 


Suppose now that the assumption (77) is dropped, and consider the special case 
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For this value, the square root in equation (76) is so small that the roots are 
coincident. These roots are both 


h ig 
and the corresponding modal ratios (74) are both 


showing that the amplitudes of the motion in the two principal modes are equal and 
in quadrature; the free motion of the particle therefore becomes a circle. The 
conjugate value of A, changes the sign of the result (84), so that a second circular 
motion, in the opposite direction to the first, is also possible. The frequency for 
these motions is the imaginary part of A, namely 


which is midway between the two natural frequencies which exist when b,, is zero. 


A further special case can be examined when 


a,o, 
This gives the roots of (76) as 
h ie b,, 
0, 4¢, * 
Taking the positive sign, the modal shape is 
Il, 
il, =-|], ‘ (88) 
g 
While the frequency is still iw, (1 + ia) (89) 


the damping component of A, is reduced from - b,,/2a, to 


1 
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With the negative sign in equation (87), the modal shape becomes 


Il, 

+1, ‘ ‘ (90) 
so that there are now two linear motions which are perpendicular. The frequency 
for the second motion is identical with that of the first, but the damping component 
is raised from — b,,/2a, to 


- (Pte). 


To sum up these results, using the particle system for illustration, suppose first 
that b,, is zero, so that the particle has two linear perpendicular motions whose 


frequencies are in the ratio 
(14 =) 1. 


Introduction of a small coefficient b,, draws the frequencies together slightly and 
transforms the motions into thin ellipses whose axes coincide with the original 
perpendicular motions. If b,, is increased further, the frequencies are drawn closer 
until eventually they coincide. The motion has by then been transformed into a 
circle and motion in either direction round it is possible. If 5,, is increased still 
further, the motions become linear again, but they are at 45° to the original linear 
motions. The damping of one of them is increased and of the other is reduced. 


All these results will be modified if the coefficients b,, and b,, differ, and in 
any particular case the motion can be found by using the preceding methods. 
In this case elliptical motions whose axes are neither coincident with, nor equally 
inclined to, the principal mode directions are possible. 
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REVIEW 


Review 


“The Numerical Solution of Two-Point Boundary Problems in Ordinary Differential 
Equations,” L. Fox. Clarendon Press, Oxford, 1957. 371 pp. 60s. 


This book is one of the Oxford University Press Monographs on Numerical Analysis 
and the production is of the excellent quality which is now taken for granted from 
that publisher. It is a large book, of some 370 pages, and, before examining it, one 
may be forgiven for wondering how a topic so apparently restricted could maintain a 
reader’s interest at such length. The explanation is that the author has included detailed 
discussion of several other topics, knowledge of which is a pre-requisite for the computor 
who wishes to use the main methods of the book in order to find solutions to ordinary 
differential equations. Thus, of its ten chapters, the first is devoted to the standard 
theories of finite-differences and the second to the solution of simultaneous algebraic 
equations, while at a later stage certain properties of matrices and the theory of 
Rayleigh’s Principle are discussed with care. 


The author makes a pronounced distinction between boundary-value techniques 
and initial-value techniques, which may be broadly distinguished by saying that the 
former find values simultaneously of a dependent function at a number of points 
distributed over a range of integration, whereas the latter find values one at a time in 
succession. Problems in differential equations may also be classified as of boundary- 
value or initial-value type, being of the second type only when all the associated end- 
conditions are specified at one and the same point of the range. Not the least interesting 
aspect of the book is the author’s demonstration of how both techniques of solution 
can be adapted to each type of problem. The basic techniques are shown to admit of 
many interesting variations: one might perhaps have wished that the author could 
have assessed their relative advantages a little more definitely, for he has clearly an 
experience wider than most of his readers can ever hope to acquire. 


Running through the book, which deals in the main with first-, second-, third- and 
fourth-order ordinary equations and with eigenvalue problems governed by such 
equations, is an emphasis on the proper use of the “ finite-difference correction ”’—a 
technique which the author is himself largely responsible for developing. His description 
of this technique is eminently clear and its advantages in several applications are well 
exemplified. To become proficient in these, as in all, numerical techniques, the student 
must gain his own practical experience, and from this viewpoint it must be regretted 
that the author does not provide some simple examples, at the end of each chapter, 
for his readers to work for themselves. 


D. N. DE G. ALLEN. 
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